高云峰

(清華大學航天航空學院,北京100084)

一般來說，在理論力學的動力學教學中，絕大部分問題都只需要學生會受力分析、列寫動力學方程，并不要求進一步求出具體的運動[1]。原因之一是絕大部分問題沒有解析解，用傳統(tǒng)數(shù)學方法不好解；原因之二是學時限制，沒有更多時間教學生如何求解。

隨著計算軟件的發(fā)展，目前已經(jīng)可以很輕松地求出動力學方程的解，包括運動和力隨時間的變化關(guān)系，也可以很容易在屏幕上進行動畫演示。

本篇先從動力學中最簡單的單擺問題開始介紹：單擺的動力學方程如何求解？數(shù)值解的精度如何？大角度情況下周期如何變化？然后介紹橢圓擺的相關(guān)問題。

1 單擺的相關(guān)問題[2]

案例1：假設(shè)單擺中小球A質(zhì)量為m，尺寸不計，繩子OA長度為l，不計質(zhì)量(圖1)。試比較不同初始角度對運動的影響，以及如何證明計算的結(jié)果可靠？

圖1 單擺模型

單擺運動時，根據(jù)其受力圖(圖2)，可以在切向方向列寫系統(tǒng)的動力學方程，為

圖2 單擺的受力圖

雖然方程(1)在線性化時有解析解，但是大角度時比較復(fù)雜。但在Matlab中，可以直接調(diào)用ode45函數(shù)統(tǒng)一求解這類常微分方程，其格式是[3]

其中黑體是固定的格式，斜體是變量或自己命名的函數(shù)。其中options是積分的選項，可以控制積分的精度，′RelTol′表示積分中的相對誤差，′AbsTol′是積分中的絕對誤差。在Matlab數(shù)值積分計算中，通常number1和number2可以選為10?8～10?10(精度太高會花費更多計算時間，也沒有必要),如果省略則默認為10?3。ode45函數(shù)則是求解常微分方程的一種常用函數(shù)，其中[t,y]是積分的時間和結(jié)果；動力學方程在子程序filename中描述；y0是初始條件；[starttime:steptime:endtime]表示積分時按等步長steptime從開始積分到結(jié)束(等步長積分是為了動畫演示時每一幀時長相同)。

注意ode45求解標準的一階常微分方程組，因此要把動力學中的二階微分方程轉(zhuǎn)換為一階微分方程組。把方程(1)這樣處理：設(shè)y1=θ,y2=，得到一階微分方程組

初始條件為y1(0)=θ(0)，y2(0)=˙θ(0)。

單擺動力學方程求解的源代碼見圖3，子程序源代碼見圖4。

圖3 單擺動力學求解的源代碼

圖4 動力學子程序源代碼

主程序源代碼中g(shù)lobal是表示全局變量，在主程序中賦值后在其他子程序中可以直接調(diào)用；plot是畫線段的函數(shù)，如果很密集，各段微小的線段就構(gòu)成了曲線；y(:,1)表示y數(shù)組中的第1列，實際上就是每一步計算得到的單擺角度；xlabel和ylabel是給圖中坐標軸加上標注，科學論文中一般需要知道坐標是什么含義，什么單位。

子程序用function開頭，注意在Matlab中子程序文件名與函數(shù)名相同；zeros(2,1)表示生成一個2×1的列陣，里面元素都是0；子程序中的動力學方程直接按照式(3)寫出。

數(shù)值計算中初始角度可以變化，其他參數(shù)不變，如下所示

圖5是方程(1)在不同初始角度下的解(暫時沒有考慮空氣阻力)，可以看到解的周期與初始參數(shù)有關(guān)，這也是非線性方程的特點。具體計算時可以把不同的初值積分結(jié)果賦值給y1，y2，y3，y4，在畫圖plot命令后使用hold on命令可以把不同的曲線疊加在一起；title是給圖命名；legend相當于圖例，可以給不同的曲線命名以示區(qū)別，見圖6。

圖5 不同角度單擺的解

圖6 多個曲線畫在同一圖上

2 數(shù)值計算的可靠性

數(shù)值計算通常是有誤差的，包括數(shù)值的截斷誤差、算法本身的誤差等。不過隨著計算技術(shù)的發(fā)展，可顯示的最小值越來越小，例如在Matlab中輸入realmin(′single′)，顯示出單精度最小正浮點數(shù)為1.175 494 4×10?38，所以普通計算精度完全滿足要求。

那么方程(1)的數(shù)值積分精度如何呢？這首先取決于options中的精度控制，在前面所說的10?8～10?10情況下，積分計算精度也基本是這一數(shù)量級。怎么證明這一點呢？可以通過方程的守恒量來判斷。

不考慮空氣阻力時，方程(1)有守恒量(機械能)

守恒量理論上是一個常數(shù)，因此可以用于檢測數(shù)值計算的可靠性。在數(shù)值計算中，可以先把每一步的角度、角速度等量計算出來，然后計算每一步的E?？紤]到計算有誤差，這個“守恒量”應(yīng)該在極小范圍內(nèi)變化才合理。圖7是不同初始條件下方程的積分結(jié)果，看起來很平坦，從而證明積分的結(jié)果很可靠。

圖7 不同條件下的積分常數(shù)

如果想了解積分常數(shù)變化的細節(jié)，利用能量之差是有效的方法。能量之差定義為

如果沒有事先指定等比例(用axis equal命令表示x和y軸等比例)，Matlab在畫圖時坐標軸會自動根據(jù)數(shù)據(jù)范圍進行縮放，因此能量之差的細節(jié)變化就可以反映出來(圖8)，根據(jù)能量之差，可以看出系統(tǒng)機械能的變化范圍與積分的允許誤差范圍同階，初始角度小時計算誤差更小。

圖8 能量之差

上述分析表明數(shù)值計算是可靠的，然后再考慮系統(tǒng)有阻尼的情況。圖9顯示了不同阻尼情況下擺角的變化情況，初始擺角均是30°，其中β=n/(2m)是阻尼系數(shù)。圖10顯示了阻尼對系統(tǒng)機械能的影響，最后能量趨于平坦，為系統(tǒng)靜止時的機械能。圖10中一個細節(jié)是：阻尼系數(shù)較小時，機械能呈現(xiàn)臺階狀的下降，這是因為擺角每次到最大幅值時速度接近，此時空氣阻力很小，機械能接近守恒，所以機械能在下降過程中就形成了臺階。

圖9 阻尼對擺角的影響

圖10 阻尼對機械能的影響

3 橢圓擺的相關(guān)問題

橢圓擺是動力學中的一個典型問題，我們想了解一下：橢圓擺的運動是周期的嗎？其周期是多少？

案例2：橢圓擺(圖11)由質(zhì)量為mA的滑塊A和質(zhì)量為mB的單擺B構(gòu)成?；瑝K可沿光滑水平面滑動，AB桿長為l，質(zhì)量不計。試建立系統(tǒng)的運動微分方程，并求水平面對滑塊A的約束力。

圖11 橢圓擺模型

系統(tǒng)有2個自由度，建立Oxy坐標系，取A點位移x和AB桿的相對轉(zhuǎn)角φ為廣義坐標(圖12)，傳統(tǒng)分析可得到系統(tǒng)的運動微分方程以及壓力的表達式(過程略)為

通常理論力學教材或教學到了這一步就算結(jié)束了。橢圓擺在運動過程中A滑塊和AB桿是否為周期運動？是否為正弦運動？壓力是如何變化的？類似這樣的問題傳統(tǒng)方法都不好回答。

系統(tǒng)的運動是周期的嗎？在式(7)中消去¨x，有

可以看出，擺角的方程(9)在小角度、線性化后是周期的；類比單擺，大角度情況下擺動周期與初始條件有關(guān)。下面的計算只改變橢圓擺的初始角度φ0，其他參數(shù)均不變。

圖13顯示了初始條件對擺角的影響，看起來不同條件下擺角都是周期函數(shù)，但是周期不同。圖14進一步顯示了周期和初始條件的關(guān)系，同時表明：初始角為小量時與線性化的結(jié)果比較接近。

圖13 橢圓擺的擺角曲線

圖14 橢圓擺的擺動周期

圖15和圖16表明位移和壓力也是周期函數(shù)，但是初始大角度時位移和壓力曲線已經(jīng)明顯偏離標準正弦曲線了。如果沒有數(shù)值計算，這些細節(jié)無法得到。

圖15 橢圓擺的位移曲線

圖16 橢圓擺的壓力曲線

另外還有一個現(xiàn)象，系統(tǒng)的位移、擺角、壓力雖然都是周期函數(shù)，但是周期并不相同。雖然橢圓擺的擺角看上去像余弦曲線(圖13)，但是否為標準的余弦曲線？這可以使用快速傅里葉變換來分析。在Matlab中，可以直接調(diào)用fft(dataname)函數(shù)求出數(shù)據(jù)dataname的頻率。

在處理橢圓擺之前，先看看fft的效果，設(shè)測試函數(shù)為

其中A0=0.5,A1=1,A2=0.4,A3=0.2;f1=3,f2=5,f3=10.6，其時域圖(圖17)看不出什么規(guī)律，但經(jīng)過快速傅里葉變換變化后，在頻域圖中就可以清楚看出只會在f=fi處有一個孤立的峰值A(chǔ)i，而在其他位置均為0：因此在頻域圖中可以直接得出測試函數(shù)(11)中的所有參數(shù)(圖18)。

圖17 測試函數(shù)的時域曲線

圖18 測試函數(shù)的頻域曲線

下面對橢圓擺的擺角和位移進行快速傅里葉變換變化(位移通過平移去掉直流分量)，為了方便比較不同初始條件的結(jié)果，把快速傅里葉變換變化后的結(jié)果歸一化處理：最大值取為1，處理后的結(jié)果見圖19和圖20，可以發(fā)現(xiàn)：

圖19 角度曲線的頻譜

圖20 位移曲線的頻譜

(1)可以看到各條曲線都有一個明顯的峰值，但其他位置還有連續(xù)的不全為0的值，且大角度時高倍頻處還會有小峰值，這說明橢圓擺的擺角不再是標準的余弦函數(shù)；

(2)由于峰值相對明顯，所以在時域圖中看起來很像余弦曲線；

(3)各曲線最大峰值對應(yīng)的主頻率不同：頻率隨初始角度增加而減少，或周期隨初始角度增加而增加，符合圖14的結(jié)論，這反映了非線性函數(shù)的周期或頻率與初始條件有關(guān)。

4 小結(jié)

現(xiàn)代科學的發(fā)展，使得計算的重要性大為提升，它和理論分析、試驗驗證同為科學研究的三大重要手段。例如，1963年MIT的氣象學家Lorenz在計算大氣對流問題時發(fā)現(xiàn)了混沌現(xiàn)象、現(xiàn)代飛機的設(shè)計涉及大量動力學軟件的計算，等等。

在這種背景下，在動力學中引入Matlab，除了會列寫動力學方程，還能計算和演示系統(tǒng)整個運動過程，便于發(fā)現(xiàn)系統(tǒng)豐富的動力學現(xiàn)象。

本篇著重介紹了動力學中運動微分方程的求解，首先把二階微分方程轉(zhuǎn)換為一階微分方程組，然后采用榮格庫塔法進行求解。涉及到的Matlab核心函數(shù)是odeset(設(shè)置積分精度)和ode45(求微分方程)；plot和hold on命令可以實現(xiàn)多個圖疊在一起。另外介紹了fft函數(shù)(把時域數(shù)據(jù)轉(zhuǎn)換到頻域)，可以分析復(fù)雜曲線的頻率或周期。

可以看出，借助Matlab可以更加深入地研究動力學問題，解決傳統(tǒng)分析方法無法處理的問題。