基于MATLAB單擺運(yùn)動(dòng)的數(shù)值分析

馬 堃

(黃山學(xué)院信息工程學(xué)院，安徽黃山245041)

單擺作為簡(jiǎn)諧振動(dòng)的一種物理模型，幾乎在每一本大學(xué)物理教材中都會(huì)介紹[1-2]。但由于其動(dòng)力學(xué)方程的非線性特點(diǎn)，在教學(xué)過(guò)程中往往只研究其小擺角的情況，此時(shí)動(dòng)力學(xué)方程簡(jiǎn)化為線性方程，其解為一余弦函數(shù)。而對(duì)于大擺角單擺的運(yùn)動(dòng)狀態(tài)隨時(shí)間的演化規(guī)律，如何在大學(xué)本科階段討論單擺振動(dòng)時(shí)的物理規(guī)律是大學(xué)物理力學(xué)部分一個(gè)教學(xué)難點(diǎn)。早在1984年，趙炳林先生就對(duì)單擺振動(dòng)的嚴(yán)格周期進(jìn)行了分析，將積分方程中的被積函數(shù)進(jìn)行級(jí)數(shù)展開(kāi)，從而得到精確的級(jí)數(shù)解[3]；隨后，人們探索了多種近似算法，如線性插值[4]、格林函數(shù)[5]、沖擊波解[6]等。隨著計(jì)算機(jī)的普及以及數(shù)學(xué)計(jì)算軟件的發(fā)展，人們?cè)诮虒W(xué)過(guò)程中逐漸探索利用數(shù)學(xué)計(jì)算軟件輔助復(fù)雜公式和方程的求解[7-8]。然而，這些研究主要集中在單擺的周期計(jì)算上，對(duì)非線性單擺的擺動(dòng)和能量轉(zhuǎn)化等情況討論的較少。本文利用MATLAB軟件對(duì)單擺非線性運(yùn)動(dòng)進(jìn)行研究，通過(guò)線性近似處理，討論非線性單擺與線性單擺運(yùn)動(dòng)狀態(tài)、能量以及周期隨時(shí)間演化的特點(diǎn)。

1 理論模型

設(shè)一單擺的擺長(zhǎng)為l，擺球質(zhì)量為m，阻尼系數(shù)為γ，則其滿(mǎn)足的動(dòng)力學(xué)微分方程可以寫(xiě)成為

其中θ是單擺偏離平衡位置的角位移，上式是一個(gè)非線性微分方程，該方程無(wú)嚴(yán)格的解析解。當(dāng)單擺做小角度擺動(dòng)時(shí)，有sinθ≈θ關(guān)系，從而方程(1)可以簡(jiǎn)化為一個(gè)線性微分方程，即

當(dāng)γ=0時(shí)，即無(wú)阻尼振動(dòng)，(1)式進(jìn)一步退化為

上式即為線性單擺所滿(mǎn)足的動(dòng)力學(xué)方程，該方程及其解在一般的教科書(shū)上都有。由于微分方程(1)式中非線性項(xiàng)sinθ的存在，無(wú)法給出其解析解。為此，本文采用數(shù)值求解方法中的四階龍格-庫(kù)塔法[9]對(duì)方程進(jìn)行求解。首先對(duì)二階微分方程作降階處理，令?，則(1)式可以降階為一階微分方程組，即

對(duì)于上式，在給定初始條件y1(t=0)=θ0，y2(t=0)=0的情況下，可以按照四階龍格-庫(kù)塔法的計(jì)算方法，讀者可以自行編寫(xiě)計(jì)算程序，也可以直接調(diào)用現(xiàn)有的程序包進(jìn)行計(jì)算，本文以MATLAB軟件為例，調(diào)用“ode45”命令進(jìn)行求解，具體編寫(xiě)程序如下：

F=@(t，y)[y(2);-g/L*y(1)-gamma/2/m*y(2)]%定義方程組(4)式

[t，y]=ode45(F，[0：0.001：20]，[theta0，v0]);%利用四階龍格-庫(kù)塔法求解，解賦值給[t，y]

本文將利用上面的MATLAB程序分別對(duì)線性擺和非線性擺進(jìn)行數(shù)值計(jì)算，并對(duì)計(jì)算結(jié)果進(jìn)行分析。方便起見(jiàn)，在計(jì)算中取m=1，l=1，g=9.8，均為國(guó)際單位。

2 線性近似下的單擺振動(dòng)

2.1 運(yùn)動(dòng)學(xué)方程

在線性近似下，只需要將一階微分方程組(4)中的siny1用y1代替即可。初始角位移和初始角速度分別取θ0=π/36和ω0=0，利用四階龍格-庫(kù)塔法可以得到單擺線性振動(dòng)時(shí)的數(shù)值解，分別是角位移y1=θ和角速度。圖1給出了阻尼系數(shù)在0～1范圍內(nèi)單擺角位移隨時(shí)間的演化規(guī)律。從圖中可以看出：無(wú)阻尼時(shí)，單擺的作等振幅振動(dòng)；有阻尼時(shí)，振幅隨著時(shí)間的推移逐漸衰減，且隨著阻尼系數(shù)的增加，衰減的越來(lái)越快。同時(shí)，我們發(fā)現(xiàn)雖然在有阻尼下單擺的振幅逐漸減小，但其擺動(dòng)的周期仍然為一定值。

圖1 線性近似下單擺作擺動(dòng)的角位移

2.2 機(jī)械能

根據(jù)角位移y1=θ和角速度y2=θ?，可以進(jìn)一步得出單擺擺動(dòng)過(guò)程中能量隨時(shí)間的演化情況。取單擺擺動(dòng)最低點(diǎn)時(shí)為勢(shì)能零點(diǎn)，初始動(dòng)能為Ek0=0 ，則初始時(shí)刻勢(shì)能為Ep0=mgl(1-cosθ0)，初始時(shí)刻總機(jī)械能為E0=Ek0+Ep0=mgl(1-cosθ0)，在單擺擺動(dòng)過(guò)程中任意時(shí)刻的勢(shì)能為

圖2分別給出了線性單擺動(dòng)能、勢(shì)能和總機(jī)械能隨時(shí)間演化情況，其中圖2(a)是無(wú)阻尼情況，圖2(b)是阻尼γ=0.5時(shí)的情況?？梢钥闯鰟?dòng)能、勢(shì)能和機(jī)械能均在周期性變化，動(dòng)能增加，勢(shì)能減小，兩者相互轉(zhuǎn)化，周期值始終保持不變。仔細(xì)分析圖2(a)可以發(fā)現(xiàn)，系統(tǒng)的總機(jī)械能會(huì)隨著時(shí)間演化而稍有波動(dòng)，且隨著擺角的增加波動(dòng)的越大?？倷C(jī)械能的波動(dòng)意味著機(jī)械能不再“守恒”。當(dāng)然，這并不是真的機(jī)械能不守恒，而是由于對(duì)微分方程(1)作了線性近似導(dǎo)致的。

圖2 線性近似下單擺的機(jī)械能

3 非線性單擺的振動(dòng)

3.1 運(yùn)動(dòng)學(xué)方程

下面我們將采用嚴(yán)格的數(shù)值求解，給出非線性微分方程(1)式的解，并討論非線性單擺運(yùn)動(dòng)過(guò)程中角位移和機(jī)械能等物理量的演化規(guī)律。

圖3 線性近似和非線性擺動(dòng)時(shí)角位移對(duì)比

圖3給出了初始角位移θ0=π/3無(wú)阻尼單擺振動(dòng)的角位移隨時(shí)間的變化關(guān)系，其中實(shí)線是非線性求解的結(jié)果，點(diǎn)劃線是線性近似下的結(jié)果?？梢钥闯?，無(wú)阻尼時(shí)，無(wú)論是線性單擺，還是非線性單擺均做等振幅擺動(dòng)，非線性擺動(dòng)周期比線性擺動(dòng)的周期大，即非線性擺隨著時(shí)間的推移擺動(dòng)越來(lái)越慢。

圖4 非線性單擺的角位移

圖4給出了不同初始角位移下非線性單擺的角位移隨時(shí)間演化情況，其中圖4(a)是無(wú)阻尼擺動(dòng)，圖4(b)是阻尼擺動(dòng)，其阻尼系數(shù)γ=0.5。從圖中可以看出，振動(dòng)周期隨著初始角位移的增加而增大，周期隨著時(shí)間推移逐漸減小，且在相同初始角位移情況下，無(wú)阻尼振動(dòng)的周期大于有阻尼振動(dòng)的周期。這是由于系統(tǒng)受到外界阻尼的影響，能量逐漸耗散，單擺擺動(dòng)的振幅也越來(lái)越小，阻尼系數(shù)越大，振幅衰減的越快。

3.2 機(jī)械能

圖5給出非線性單擺動(dòng)能、勢(shì)能和機(jī)械能隨時(shí)間演化的情況，其中圖5(a)是無(wú)阻尼擺動(dòng)，圖5(b)是阻尼擺動(dòng)。無(wú)阻尼時(shí)，動(dòng)能和勢(shì)能相互轉(zhuǎn)化，此消彼長(zhǎng)，但系統(tǒng)的總機(jī)械能保持不變，等于系統(tǒng)初始的機(jī)械能，滿(mǎn)足能量守恒定律，這一點(diǎn)與線性擺不同。有阻尼時(shí)，動(dòng)能和勢(shì)能與會(huì)隨著時(shí)間相互轉(zhuǎn)換，但由于阻尼的存在，在轉(zhuǎn)換的過(guò)程中，總機(jī)械能隨著時(shí)間推移而逐漸較小。

圖5 非線性單擺的機(jī)械能

3.3 周期

可以看出，線性單擺的周期只與擺長(zhǎng)有關(guān)，與初始角位移無(wú)關(guān)，而非線性單擺的周期不僅與擺長(zhǎng)有關(guān)，還與初始的角位移有關(guān)。表1給出了不同初始角位移下非線性單擺的周期，其中T表示本文數(shù)值計(jì)算得到的結(jié)果，T0表示線性擺周期，Ref[T/T0]表示(8)式計(jì)算的結(jié)果。從該表中可以看出，隨著初始角位移的增大，單擺的周期增大，且阻尼振動(dòng)的周期小于無(wú)阻尼時(shí)的周期。通過(guò)與嚴(yán)格的積分形式解(8)式的結(jié)果比較可見(jiàn)，本文數(shù)值計(jì)算的結(jié)果與(8)式結(jié)果的誤差均在0.3%之內(nèi)。

表1 不同初始角位移下非線性擺的周期

圖6給出了單擺的擺動(dòng)周期與初始角位移之間的關(guān)系，其中虛線表示線性近似下的結(jié)果。在線性近似下，單擺的周期與初始角位移無(wú)關(guān)，是一個(gè)定值。圖中可以看出非線性擺的周期與不僅與初始角位移有關(guān)，而且與阻尼有關(guān)。初始角位移越大，周期越大，阻尼越大，周期越小。

圖6 單擺周期與初始角位移的關(guān)系

4 總結(jié)

本文利用四階龍格-庫(kù)塔法嚴(yán)格求解單擺的非線性動(dòng)力學(xué)微分方程，討論了非線性單擺的角位移、能量和周期隨時(shí)間的演化規(guī)律。在非線性解下，單擺的機(jī)械能不再隨著時(shí)間發(fā)生變化。雖然本文的分析是基于數(shù)值求解方法，但是隨著計(jì)算機(jī)的普及以及數(shù)學(xué)計(jì)算軟件的發(fā)展，我們認(rèn)為有必要、也有條件在大學(xué)本科階段普及數(shù)值計(jì)算方法，以幫助加深對(duì)物理概念的理解，拓展學(xué)生運(yùn)用所學(xué)知識(shí)解決實(shí)際問(wèn)題的能力。