單擺問題の近似解

01-26

有關單擺的知識，我們早在中學階段就已經接觸過了。但是，你可曾想過，如果單擺的上懸點運動起來，那會怎樣(⊙_⊙)？今天我們就一起來看一下這個有趣の問題吧(*^__^*) ......

首先，需要說明的是：

1，空氣阻力雖然很小很小，但依然存在，是絕對不可以忽略的，我們假設阻力大小和速度大小成正比。

2，我們假設重物要比繩子重得多得多，而且單擺運動時繩子始終緊繃（當然這個假設在某些情況下是不成立的，我們在這裡僅僅討論假設成立情況下單擺的運動）。

第一部分：上懸點固定の情況

1，以上懸點為原點，水平方向為x軸，豎直方向為z軸，建立平面直角坐標系。很顯然，既然繩子一直緊繃，那麼我們就可以用一個約束方程來表徵繩子對重物運動的影響，約束方程為： $varphi left( x,z ight) =x^{2} +z^{2} -l^{2} =0$ ，其中 $l$ 為繩子的長度。

2，下面我們要建立描述單擺運動的微分方程，體系的拉格朗日函數為： $L=frac{1}{2} mleft( x^{2} +z^{2} ight) -mgz$ ，但現在繩子の約束也會影響重物的運動。可以看出，我們要處理的是一個含有約束條件的力學問題，可以利用拉格朗日乘數法構造輔助泛函進行求解。輔助泛函為： $ilde{L} =L+lambda left( t ight) varphi$ ，其中 $lambda left( t ight)$ 是拉格朗日乘子。下面要做的就是把輔助泛函帶入拉格朗日方程，就可以得到體系的微分方程了。

3，將 $ilde{L} =L+lambda left( t ight) varphi$ l帶入拉格朗日方程並化簡得： $x+alpha ^{2} x=0$ ， $z+alpha ^{2} z=-g$ ，其中 $alpha ^{2} =-frac{2lambda left( t ight) }{m}$ 。從數學上看， $alpha$ 是拉格朗日乘子 $lambda left( t ight)$ 變形而來，從數學上分析它是十分困難的，所以我們從物理角度對他進行分析。對 $alpha$ 做量綱分析並結合我們之前的知識就可以發現， $alpha$ 就是是單擺的固有角頻率。

4，對於單擺而言，它的固有角頻率 $alpha$ 是取決於擺角 $heta$ 的，具體表達式為： $alpha =frac{2pi }{T}$ ， $T=4left( frac{l}{g} ight) ^{frac{1}{2} } int_{0}^{frac{pi }{2} } frac{1}{sqrt{1-sin^{2} frac{ heta }{2}sin^{2} varphi } } dvarphi$ ，我們對它進行冪級數展開，取常數項得： $alpha =left( frac{g}{l} ight) ^{frac{1}{2} }$ ，這正是我們做近似處理的地方。當擺角 $heta =10^{circ }$ 時，相對誤差為0.154%；當擺角 $heta =20^{circ }$ 時，相對誤差為0.689%；當擺角 $heta =30^{circ }$ 時，相對誤差為1.599%。可見，擺角不大時，這種近似的近似程度是非常好的。

5，加入線性的空氣阻力項後，體系的微分方程為：

$x+hx+alpha ^{2} x=0$ ， $z+hz+alpha ^{2}z=-g$ ，方程的解為：

$x(t)=A_{0} e^{-eta t} cos(omega _{f} t+varphi _{0})$ ， $z(t)=C_{0} e^{-eta t} cos(omega _{f} t+varphi _{0}) -frac{g}{alpha ^{2} }$ ，其中， $eta =frac{h}{2} ,omega _{f} =sqrt{omega _{0}^{2}- {eta }^{2} }$ ,可見，由於空氣阻力（能量耗散），重物會逐漸靜止下來。

第二部分：上懸點做任意形式的運動

1，以空間中一點為原點建立直角坐標系，設上懸點運動的參數方程為： $x_{0} =x_{0} (t);y_{0} =y_{0} (t);z_{0}=z_{0} (t)$ ，和第一部分一樣，代表繩子對重物影響的約束方程為： $varphi (x,y,z)=(x-x_{0}(t))^{2} +(y-y_{0}(t))^{2} +(z-z_{0} (t))^{2} -l^{2} =0$ ，同理用拉格朗日乘數法構造輔助泛函 $ilde{L}=L+lambda left( t ight)varphi$ ，帶入拉格朗日方程，並加入線性阻尼項得：

$x+hx+alpha ^{2} x=alpha ^{2} x_{0}(t)$ ， $y+hy+alpha ^{2} y=alpha ^{2} y_{0}(t)$ ， $z+hz+alpha ^{2} z=alpha ^{2} z_{0}(t) -g$ ，

2，其中 $alpha$ 即為單擺的固有角頻率，這個結論可以由量綱分析結合第一部分內容得到。當上懸點不動時，體系還原為第一部分中的單擺，所以， $alpha =left( frac{g}{l} ight) ^{frac{1}{2} }$ 。當然這也是一種近似。

3，可以看出，以上三個方程都是二階常係數線性非齊次の微分方程，其通解等於相應其次方程通解加上一個非齊次方程の特接，所以，只要知道上懸點運動，大部分情況下都可以求出微分方程的通解。

第三部分：上懸點做周期性運動

1，周期性運動有很多種，這裡我們只分析一種，即，上懸點在水平面內做勻速圓周運動。其他的周期性運動都可以仿照它做類似的討論。我們一上懸點軌跡的圓心為原點，水平面為xoy平面，豎直方向為z軸，建立平面直角坐標系。

2，根據第二部分的結論可知，體系的微分方程為：

$x+hx+alpha ^{2} x=alpha ^{2} Rcos(omega _{0} t)$

$y+hy+alpha ^{2} y=alpha ^{2} Rsin(omega _{0} t)$

$z+hz+alpha ^{2} z=-g$

微分方程組的通解為：

$x(t)=A_{0} e^{-eta t} cos(omega _{f} t+varphi _{0} )+Acos(omega _{0} t-varphi )$

$y(t)=B_{0}e^{-eta t} cos(omega _{f} t+varphi _{0})+Bcos(omega _{0}t-varphi )$

$z(t)=C_{0} e^{-eta t} cos(omega _{f} +varphi _{0})-frac{g}{alpha ^{2} }$

可見，開始時z方向做簡諧運動，x，y方向做准周期運動。由於空氣阻力，體系的能量不斷散失，z方向振幅最後衰減為零，x，y方向衰減為一個和上懸點角速度相同的勻速圓周運動。

3，對於方程通解中，最後「剩餘」圓周運動，其半徑為： $A=B=frac{alpha ^{2}R }{((alpha ^{2}-omega _{0}^{2} )^{2} +4eta ^{2} alpha ^{2} )^{frac{1}{2} } }$ ，其相對於上懸點的勻速圓周運動の滯后角滿足： $tanvarphi =frac{2eta omega _{0} }{alpha ^{2} -omega _{0}^{2} }$ ，

（1）由於空氣阻力，體系能量不斷散失，最後剩下の是一個勻速圓周運動（穩定軌跡）。

（2）重物做勻速圓周運動的角速度和上懸點相同。

（3）重物運動半徑與上懸點相比可能大也可能小。

（4）由於空氣阻力非常小，所以 $eta approx 0$ ，滯后角的正弦值趨近於零。當 $alpha >omega _{0}$ 時， $tanvarphi$ 正向趨近於零，滯后角 $varphi =0$ ；當 $omega _{0} >alpha$ 時， $tanvarphi$ 負向趨近於零，滯后角 $varphi =pi =180^{0}$ .

（5）當 $alpha >omega _{0}$ （滯后角為0度）時，重物運動半徑一定大於上懸點運動半徑（讓 $alpha$ 趨近於 $infty$ 即可證明）。所以，體系的穩定軌跡（勻速圓周運動）只有三種：滯後0度半徑大；滯後180度半徑小；滯後180度半徑大。

第四部分：兩種特殊情況的初步分析

在第三部分的基礎上，我們討論兩種特殊的情況，①阻力為零。②極高轉速（阻力較大）

1，阻力為零的情況：

阻力為零時，體系不是耗散系統，即，體系沒有吸引子。這種情況下我們會得到無窮多種依賴於初值的准周期運動，其龐加萊局域截面上的點組成一個閉合的曲線。

2，轉速極高的情況：

在轉速極高的情況下，方程中必會出現某種形式的非線性剛度項，和非線性阻尼項。這種情況下，我們無法給出方程的解析解。我們可以求出方程的數值解，並在相空間中繪出相空間的軌跡。這種情況下，龐加萊局域截面上會呈現出大範圍具有分形結構的密集點。此時運動是混沌的，即體系吸引子不僅僅只有參數決定，還極其敏感的依賴於體系的初值，Lyapunov指數在一定程度上可以定量描述體系對初值的敏感程度。有了以上討論，還可以進一步對體系做混沌控制等內容的討論。

第五部分：拓展

現在考慮一個問題，如果我在重物下面再「掛」一個單擺，那將會怎樣呢？如果再「掛」n個呢？

同樣還是用拉格朗日乘數法來處理這個問題，這裡直接給出n個單擺連在一起の結果。

$x_{1}+hx_{1} + alpha _{1}^{2} x_{1} +alpha _{2}^{2} x_{1} =alpha _{1}^{2} x_{0} +alpha _{2}^{2} x_{2}$ ，

$x_{2}+hx_{2} + alpha _{2}^{2} x_{2}+ alpha _{3}^{2} x_{2} =alpha _{2}^{2} x_{1}+ alpha _{3}^{2} x_{3}$ ，

.............................................

$x_{n-1} +hx_{n-1} +alpha _{n-1}^{2} x_{n-1}+ alpha _{n}^{2} x_{n-1} =alpha _{n-1}^{2} x_{n-2} +alpha _{n}^{2} x_{n}$

$x_{n}+ hx_{n}+ alpha _{n}^{2} x_{n}+ =alpha _{n}^{2} x_{n-1}$

其中，上懸點運動 $x_{0}= x_{0} (t)$ 為已知函數。同理可得y，z方向運動の微分方程組。事實上，這裡的固有角頻率依然是變數，即， $alpha =alpha(x,y,z ,t)$ 。我們同樣可以對他進行冪級數展開，並取常數項，來求得近似解。