Matlab求解理論力學(xué)問題系列(三)單擺和橢圓擺的運(yùn)動(dòng)及周期

高云峰

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

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

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

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

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

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

圖1 單擺模型

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

圖2 單擺的受力圖

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

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

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

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

單擺動(dòng)力學(xué)方程求解的源代碼見圖3，子程序源代碼見圖4。

圖3 單擺動(dòng)力學(xué)求解的源代碼

圖4 動(dòng)力學(xué)子程序源代碼

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

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

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

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

圖5 不同角度單擺的解

圖6 多個(gè)曲線畫在同一圖上

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

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

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

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

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

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

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

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

圖8 能量之差

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

圖9 阻尼對(duì)擺角的影響

圖10 阻尼對(duì)機(jī)械能的影響

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

橢圓擺是動(dòng)力學(xué)中的一個(gè)典型問題，我們想了解一下：橢圓擺的運(yùn)動(dòng)是周期的嗎？其周期是多少？

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

圖11 橢圓擺模型

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

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

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

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

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

圖13 橢圓擺的擺角曲線

圖14 橢圓擺的擺動(dòng)周期

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

圖15 橢圓擺的位移曲線

圖16 橢圓擺的壓力曲線

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

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

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

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

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

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

圖19 角度曲線的頻譜

圖20 位移曲線的頻譜

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

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

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

4 小結(jié)

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

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

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

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