（3）根據數學模型建立計算機仿真模型（編程） 計算機仿真就是對數學模型的數值求解。我們將微分方程進行形式上的 變換以便於數值 求解。由（1.2）和（1.3）式得 v(t+ dt) =v(t) + dv=v(t) +adt (1.9) 和 s(t+ dt) =s(t) + ds=s(t) +v(t)dt (1.10) 注意，這種變形僅僅是將方程轉換為一種在自變量（時間）上的「遞 推」表達式，並沒有進行 解析求解。利用（1.9）和（1.10），在已知當前時刻t 的瞬時位移，瞬 時速度和加速度的情況下，我們就可以推知下一個無限鄰近的時刻t+ dt 上物體新的瞬時位移，瞬時速度和加速度，這也就是微分方程數值求解 的基本思想。在數值求解中，無窮小量dt 需要用一個很小的數量¢t來近 似，¢t 稱為微分方程的數值求解步長，通常也稱為仿真步進。顯然，這 種微分方程的遞推求解總是近似的，求解精度與步長有關。下面我們就 用程序來實現這個求解過程。仿真時間範圍設置為0到2秒，為了使得仿 真計算的誤差明顯一些，我們故意採用了較大的仿真步長¢t = 0:1。讀者 可以自己修改仿真步長來觀察計算精度的變化情況。我們將瞬時位移作 為仿真輸出變量，並同時通過（1.8）式計算出解析結果以相互對照。 （4）執行仿真和結果分析 g=9.8; % 重力加速度 v=0; % 設定初始速度條件 s=0; % 設定初始位移條件 t=0; % 設定起始時間 dt=0.1; % 設置計算步長 N=20; % 設置仿真遞推次數. 仿真時間等於N與dt的乘積 for k=1:N v=v+g*dt; % 計算新時刻的速度 s(k+1)=s(k)+v*dt; % 新位移 t(k+1)=t(k)+dt; % 時間更新
畫」仿真的效果。當然，這樣是以犧牲計算速度為代價的。修改後的程 序文件代碼如下。 〔程序代碼〕ch1example1prg2.m % ch1example1prg2.m g=9.8; % 重力加速度 for L=1:5 % 仿真重複5次以便於觀察 v=0; % 初始速度 s=0; % 初始位置 t=0; dt=0.01;% 計算步長 for k=1:200 v=v+g*dt; % 速度 s=s+v*dt; % 位移 t=t+dt; % 時間 plot(0,-s,'o'); axis([-2 2 -20 0]); % 坐標範圍固定 text(0.5,-1,['当前时间：t=',num2str(t)]); text(0.5,-2,['当前速度：v=',num2str(v)]); text(0.5,-3,['当前位置：s=',num2str(s)]); set(gcf,'DoubleBuffer','on');% 雙緩衝避免作圖閃爍 drawnow; % 立即作圖 end end 在程序中，計算步進重新設置為0:01，並且重複仿真多次以便於演示。 在作圖語句之後設 置了固定的顯示坐標範圍，並通過text語句在圖上動態顯示當前計算時 刻的結果值。語句 set(gcf,'DoubleBuffer','on'); 將顯示設備雙緩衝開啟，以避免屏幕刷新引起的閃爍。而 drawnow 是立即執行繪圖指令。在仿真中途可按鍵Ctrl+C來終止程序執行。程序 執行過程中將顯示出物體墜落的動畫效果，其中的一個幀如圖1.2所 示。
பைடு நூலகம்
用Matlab常微分方程求解器重新仿真實例1.2的乒乓球彈跳模型，要求也 能 夠在仿真過程中動態顯示球的墜落和反彈過程，在仿真結束後輸出小球 位移和速度的變化曲 線。設小球的初始位置距離水平面1米，由靜止狀態開始自由墜落，並 設反彈瞬間的速度衰減 係數為K= 0:85。 數學模型與實例1.2完全相同，為了利用標準求解器進行求解，首先將 模型中的方程改寫 為標準的狀態方程形式，小球在空中時，有運動方程
並注意小球碰撞水平面瞬間（在反彈之前）的條件是y60且v60。碰撞瞬 間速度發生反向 並以係數K衰減。據此編寫程序，程序中以x1(t)表示速度狀態變量v(t)， 以x2(t)表示位移 狀態變量y(t)，並以矩陣形式描述，代碼如下。 v0=0; y0=1; %球的初始狀態 x_state=[v0,y0];%將初始狀態賦值到狀態變量中
下面我們通過對自由落體的仿真試驗來說明計算機仿真的過程。 〔實例1.1〕試對空氣中在重力作用下不同質量物體的下落過程進行建 模和仿真。已知重 力加速度g= 9:8m/s^2，在初始時刻t0= 0s時物體由靜止開始墜落。空氣 對落體的影響可以 忽略不計。 （1）建立數學模型 首先，我們根據物理原理建立自由落體的數學模型。在空氣阻力可忽略 不計時，質量為m 的物體在自由墜落過程中受到豎直向下的恆定重力的作用，由牛頓第二 定律，我們知道，重 力F，加速度a以及物體質量m之間的關係是： F=ma (1.1) 其中加速度就是重力加速度，即a=g。根據題設，初始時刻為t0= 0，物 體的初始速度 為v(t0) = 0，並設物體下落的瞬時速度為v(t)。設物體在t 時刻的位移為 s(t)，並設初始位移為零，即s(t0) = 0。根據加速度、速度、位移三者之 間的微積分關係，我們得到一組數學方程 a =dv/dt v =ds/dt F = ma 以及初始條件（也稱為方程的邊界條件） v(t0) = 0 (1.5) s(t0) = 0 m (1.6 我們需要得出不同時刻物體的運動狀態，即物體的瞬時速度和瞬時位 移。至此，我們得到 了自由落體的數學描述。 （2）數學模型的解析分析 數學模型建立之後，可以嘗試對其進行解析求解。解析結果可以幫助我 們驗證仿真數值 結果。對於這個數學模型，其求解十分簡單。只要對加速度方程、速度 方程進行積分並代入初始條件，就能得到我們熟知的結果：
實例1.2〕對乒乓球的彈跳過程進行仿真。忽略空氣對球的影響，乒乓 球垂直下落，落點 為光滑的水平面，乒乓球接觸落點立即反彈。如果不考慮彈跳中的能量 損耗，則反彈前後的瞬時速率不變，但方向相反。如果考慮撞擊損耗， 則反彈速率有所降低。我們希望通過仿真得出乒乓球位移隨時間變化的 關係曲線，並進行彈跳過程的「實時」動畫顯示。 （1）數學模型 首先對乒乓球彈跳過程進行一些理想化假設。設球是剛性的，質量為 m，垂直下落。碰 擊面為水平光滑平面。在理想情況下碰擊無能量損耗。如果考慮碰擊面 損耗，則碰擊前後速 度方向相反，大小按比例係數K;0< K<1下降。在t 時刻的速度設為 v=v(t)，位移設 為y=y(t)，並以碰擊點為坐標原點，水平方向為坐標橫軸建立直角坐標 系。球體的速度以豎 直向上方向為正方向。重力加速度為g= 9:8m=s2。初始條件假設：設初 始時刻t0= 0球體的初始速度為v0=v(t0)，初始位移為y0=y(t0)。受力分 析：在空中時小球受重力F=mg作用，其中，g=-dv/dt。則在t+ dt 時刻小 球的速度為（注意其中負號是考慮了速度的方向） v(t+ dt) =v(t)¡gdt (1.11) 在t+ dt 時刻小球的位移為 y(t+ dt) =y(t) +v(t)dt (1.12) 在小球撞擊水平面的瞬間，即y(t) = 0的時刻，它的速度方向改變，大小 按比例K衰減。 當K= 1時，就是無損耗彈跳情況。因此，小球反彈瞬間（t+ dt 時刻）的
function xdot=ch2example6statefun(t, x, flag) % 乒乓球彈跳模型的標準狀態方程 % x(1)為小球速度，x(2)為小球位移 xdot=zeros(2,1); % 狀態變量矩陣初始化 xdot(1)=-9.8; % 速度加速度方程 xdot(2)=x(1); % 位移速度方程
end % 理論計算, 以便與仿真結果對照 t_theory=0:0.01:N*dt; % 設置解析計算的時間點 v_theory=g*t_theory; % 解析計算的瞬時速度 s_theory=1/2*g*t_theory.^2; % 解析計算的瞬時位移 % 作圖: 仿真結果與解析結果對比 t=0:dt:N*dt; plot(t,s,'o', t_theory,s_theory, '-'); xlabel('時間t'); ylabel('位移s'); legend('仿真結果','理論結果'); 仿真程序編寫完成並調試正確之後，運行得到的結果如圖1.1所示。從 圖中可見，仿真得 出的位移與理論結果之間存在差別，這種差別是由於微分方程數值求解 的算法和採用步長較 大而引起的。事實上，我們這裡採用的數值求解算法是最簡單的矩形積 分的方法，精度不高， 僅僅是為了說明仿真過程而已。現代仿真技術和數值計算方法中已經開 發出許多更好的微分 方程求解算法，可供直接利用。 （5）仿真程序的功能擴展 從ch1example1prg1.m的程序代碼中可見，我們採用了循環語句來實現對 微分方程的遞 推求解，每次循環就將計算時刻向前推進一個步長。全部循環執行完畢 後，就得到了一系列時刻上物體的瞬時速度和瞬時位移值，最後通過繪 圖語句將結果數據用曲線表達出來。 如果我們希望以動態方式來觀察物體墜落的過程，可以這樣設計仿真程 序：使得在數值求 解的過程中能將求解結果以圖形方式輸出出來。這樣，在數值求解不斷 更新的過程中輸出圖 形也隨之同步更新，形成一種「動畫」的效果。這種一邊計算一邊輸出 可視化結果的方式更 加形象直觀，更便於理解物理系統的工作過程，同時也方便演示，教學 講解和學術交流。 我們可以將作圖語句放在遞推計算循環內，並設置即時作圖刷新方式， 從而得到這種「動
...... legend('K=0.85','K=1'); 圖1.3中分別作出了是碰擊衰減係數K= 1和K= 0:85兩種情況下的小球彈 跳位移曲 線。對程序稍加修改就可以得到顯示小球彈跳過程的動畫。修改後的程 序文件代碼如下，讀者可運行該程序觀察不同的碰擊衰減係數下的小球 彈跳過程。 〔程序代碼〕ch1example2prg2.m % ch1example2prg2.m g=9.8; % 重力加速度 v0=0; % 初始速度 y0=1; % 初始位置 m=1; % 小球質量 t0=0; % 起始時間 K=0.85; % 彈跳的損耗係數 N=5000; % 仿真的總步進數 dt=0.005; % 仿真步長 v=v0; % 初狀態 y=y0; for k=1:N if y >0 % 小球在空中的動力方程計算 v =v -g*dt; y =y +v*dt; else % 碰擊瞬間的計算 y =-K.*v*dt; v =-K.*v-g*dt; end plot(0,y,'o'); axis([-2 2 0 1]); % 坐標範圍固定 set(gcf,'DoubleBuffer','on'); % 雙緩衝避免作圖閃爍 drawnow; end