馬 堃

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

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

1 理論模型

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

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

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

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

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

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]);%利用四階龍格-庫塔法求解，解賦值給[t，y]

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

2 線性近似下的單擺振動

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

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

圖1 線性近似下單擺作擺動的角位移

2.2 機(jī)械能

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

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

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

3 非線性單擺的振動

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

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

圖3 線性近似和非線性擺動時角位移對比

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

圖4 非線性單擺的角位移

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

3.2 機(jī)械能

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

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

3.3 周期

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

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

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

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

4 總結(jié)

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