小球在彈簧頂端的木塊上的彈性跳動問題之代碼分析

根據(jù)上次課程的原理分析及給出的代碼（http://www.yqgqt.org.cn/content/post/318856），這次課程的任務(wù)是對給出的代碼進行分析。本次課程分為三部分，一是上一節(jié)課程中的結(jié)果展示，二是代碼實現(xiàn)分析；三是代碼分析。

我們先按照給出的代碼運行程序，得出了以下圖片。

1.結(jié)果展示

圖 1.小球與彈簧塊位移圖

圖 2.小球與彈簧塊運動動畫（動圖）

得到了小球與彈簧塊的位移圖（圖1）與小球與彈簧塊運動（圖2），在圖1中看到小球與彈簧塊接觸25次，這個數(shù)值是在代碼中設(shè)定的；圖2是一張動畫圖，展示的是小球與彈簧塊的運動過程。

2.代碼實現(xiàn)分析

從上一節(jié)課程中我們了解到，小球與彈簧塊的運動過程分為兩大部分，一是未碰撞情況下，我們運用基本的受力分析以及位移，速度，加速度的關(guān)系可得方程；二是碰撞的情況下，我們采用動量守恒定律以及能量守恒定律得到方程。代碼運行應(yīng)該符合下列流程圖。

圖3 代碼實現(xiàn)流程

    首先根據(jù)給定的初值，進入方程組（2），直到小球與彈簧塊碰撞時進入方程組（4），如此便得到了小球與彈簧塊第一次碰撞間的所有過程，因為本例不存在任何阻尼的作用，所以，小球?qū)椘穑瑥椈蓧K將下降，從新進入到方程組（2），再開始一個新的循環(huán)。理論上來講，這將是一個無限的循環(huán)，但我們可以通過設(shè)置小球與彈簧塊的碰撞次數(shù)來限定，在代碼中，設(shè)定的25次。

3.代碼分析

3.1 解微分方程基本知識

    微分方程分為常微分方程（ODE）與偏微分方程（PDE），簡單來說系數(shù)為常數(shù)的為ODE，而系數(shù)為變量的為PDE。目前Matlab提供了多種解微分方程的結(jié)算器，如ode113、ode15i、ode23s、ode3t等，以后有機會將系統(tǒng)的闡述各種算法的差異以及其各種擅長的領(lǐng)域，現(xiàn)在只需要記住，一般默認ode45算法即可，Matlab公司的建議是，如果不清楚系統(tǒng)的具體情況，建議先選擇ode45。常微分解法器的輸入輸出格式有：

[t,Y] = solver(odefun,tspan,y0)；

[t,Y] = solver(odefun,tspan,y0,options)；

[t,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)；

    上述指令的“solver”所指的就是ode45、ode45, ode23, ode113, ode15s, ode23s, ode23t, or ode23tb等，其他參數(shù)的含義介紹如下：

a) odefun 用以計算微分程右邊的函數(shù)，所有解法器處理的系統(tǒng)方程形式均是y’=dydt=f(t,y)的形式，通常可由odefile模板改寫而成；

b) tspan 為積分區(qū)間的向量[t0,tf],解法器設(shè)定初值在tspan(1)，然后由tspan(1)積分至tspan(end)。如果需要得到特定時間所對應(yīng)的解，可以用tspan[t0,t1,…,tf]的形式。使用tspan向量的大小并不影響運算的時間，影響的是存儲的空間；

c) y0 為初始向量指；

d) options 為改變積分特性的參數(shù)，可以用odeset設(shè)定，如本先將積分器中的時間判定功能打開；

e) [t  Y] 為輸出值，t是時間點的行向量，Y是解答列陣，每一行Y均與時間t所對應(yīng)，每組(t,Y)所產(chǎn)生稱為事件函數(shù)，每次均要檢查是都等于0，并決定在零指時是否終止運算；

3.2   “events”函數(shù)分析

    基本形式為：[value, isterminal,direction]=events(t,y)

    其中value(i) 為函數(shù)的值，isterminal(i)=1時運算在等于0時停止，isterminal(i) =1時運算在等于0時繼續(xù)；direction(i)=1在事件函數(shù)增加時等于零，direction(i)=-1在事件函數(shù)減小時等于零。這一功能與上一小節(jié)“e）”點組合理解。

3.3   代碼分析

圖4.代碼分析1

    這一部分主要實現(xiàn)的是解微分方程，可以說是整個算法的全部。首先設(shè)定小球下落高度為50，彈簧系數(shù)為60；小球質(zhì)量為20；彈簧塊質(zhì)量為50；積分開始時間為0，結(jié)束時間為1000，注意這不是系統(tǒng)時間，1000也不是意味著1000秒，而是用于微分解算器中tspan參數(shù)設(shè)定；初值為[h0,0,0,0]，其中h0是小球初始高度，第二個“0”指的是彈簧塊的初始位置（坐標原點，為0），第三個0指的是小球的初始速度，為0意味著是靜止下降，第四個0指的是彈簧塊的初始速度。

     將tstart賦予tout，y0的轉(zhuǎn)置賦予yout，設(shè)置積分算法開啟事件判斷功能。

     接著，進入25個循環(huán)，這意味著小球與木塊將碰撞25次，下面一句語句

      [t,y,event]=ode45('xqythkfun',[tstart:0.03:tfinal],y0,options);是解算微分方程，其中涉及的'xqythkfun'利用odefile模板編寫，見下圖。

圖5 代碼分析圖2（微分方程函數(shù)）

默認計算f(t,y)，如果識別有事件發(fā)生則進入events，這也是使用子函數(shù)的形式，f（t，y）如下圖所示。

圖6 方程（2）

這就是實現(xiàn)了方程（2）的功能，events函數(shù)如下圖所示。

圖7 時間判斷函數(shù)

     這是用于判斷進入碰撞，如果碰撞了則該次積分結(jié)算器停止，判斷的條件是小球與彈簧塊接觸。direction=-1 意 味著Q值由大變到0時停止本次積分。

     接著回到圖4，為了方便查看，截取余下待分析代碼如下圖所示。

圖8 主函數(shù)中待分析代碼

其中

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

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

是為了將25次碰撞過程中的所有數(shù)據(jù)儲存，便于畫圖使用。之所以用（2：end）是因為初始條件給出了第一列的數(shù)據(jù)。

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

是為了更新y0（初始條件）狀態(tài)，由于以上已經(jīng)判斷了發(fā)生了碰撞，所以需要將碰撞后的數(shù)據(jù)作為初始條件（y0），對于位移來說，實現(xiàn)很簡單，只需要替換即可，而對于速度，需要用動量定理以及能量守恒定理來解決，就有下面的代碼：

    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);

其主要是解了方程（4）

重復(fù)25個碰撞過程，就完成了我們所設(shè)定的程序。

Tips：

1.文中所設(shè)計的方程（2）與方程（4）見上一節(jié)課程（http://www.yqgqt.org.cn/content/post/318856）；

       2.作圖將在下次課程中分析；