第一篇：Bazinga：大角度單擺的運動近似(1)：小角度近似

第二篇：Bazinga：大角度單擺的運動近似(2):周期的解與橢圓積分

上次我們已經成功解決了周期的表達式以及周期近似，接下來我們把注意力集中到微分方程 [ddot heta = - frac{g}{L}sin heta ] 的解 heta(t) 上來.此時我們不使用小角度近似，並試圖求該微分方程的數值解.在往後的幾篇文章中我會把自己的級數解法更新上來.在本文會先用一種我臨時受到啟發的解法來解決該問題.靈感出自《費曼物理講義》.

首先對該物理問題所需要的數學基礎做一個鋪墊.

圖片來自www.bing.com

本文的核心數學思想是線性近似.對於一個已知函數 y=f(x) ,由自變數在 x=a 處增加引起的增量 [Delta y = f(a + Delta x)] ,此增量可以表示為

[Delta y = f(a)Delta x + o(Delta x)]

在 [Delta x o 0] 時高階無窮小量 [o(Delta x) o 0] ，這時可以定義微分了.不過這裡我們採用近似的記法，即當 [Delta x] 很小時， [Delta y approx f(a)Delta x]

這個近似模型可以很好的運用在運動學中，比如在本文的情況下 [omega = frac{{d heta }}{{dt}}] , [alpha = frac{{domega }}{{dt}}]

所以可以近似的說，當我們的步長 [Delta t] 足夠小的時候，

[ heta ({t_n}) = heta ({t_{n - 1}}) + omega ({t_n})Delta t]

[omega ({t_n}) = omega ({t_{n - 1}}) + alpha ({t_{n - 1}})Delta t]

注意在這裡我們定義 [{t_n} = nDelta t] ,這說明了對於任意時刻的角度位移，只要知道它上一刻的角度位移和此刻的角速度就可以得到；而此刻的角速度又由上一刻的角速度和上一刻的角加速度給出.此時牛頓第二定律的微分方程給了我們便利，因為它成功的把角加速度和角度位移聯繫了起來：

[alpha ({t_n}) = frac{{{d^2}}}{{d{t^2}}} heta ({t_n}) = - frac{g}{L}sin heta ({t_n})]

這也就是說

[left{ egin{array}{l} heta ({t_n}) = heta ({t_{n - 1}}) + omega ({t_n})Delta t\ omega ({t_n}) = omega ({t_{n - 1}}) - frac{g}{L}sin heta ({t_{n - 1}})Delta t end{array} ight.]

故只要我們知道了前一時刻的角加速度和角度位移，就可以表示出這一時刻的角度位移.此時我們為了化簡計算，令常數 [frac{g}{L} = 1] .同時為了獲得較高的精確度，我們讓 [Delta t = 0.01] .

我們還需要給出初始條件.考慮到我們研究的是單擺在大角度情況下的運動，那麼我們令 t=0 時 [ heta ({t_0}) = frac{pi }{2}] ,而角速度自然是 [omega ({t_0}) = 0] .

為了畫出角度位移隨時間變化的圖像，我讓我同學 @天琴心弦 （編程大佬）幫我寫了一段Python代碼，從而幫助我畫出這個函數:

import numpy as np
import matplotlib.pyplot as plt

telement = list()
telement.append(np.pi/2)
for i in range(1000):
telement.append(0)
oelement = list()
oelement.append(0)
for i in range(1000):
oelement.append(0)

def theta(n):
if n == 0:
return telement[n]
else:
print("This is "+ str(n))
if telement[n] == 0:
telement[n] = telement[n - 1] + 0.01 * omega(n)

return telement[n]

def omega(n):
if n == 0:
return oelement[n]
else:
if oelement[n] == 0:
oelement[n] = oelement[n - 1] - np.sin(telement[n - 1]) * 0.01

return oelement[n]

x = list()
for i in range(1000):
x.append(i*0.01)
y = list()
for i in range(1000):
y.append(theta(i))

plt.plot(x, y)
plt.show()

運行之後大概是這個樣子：

我讓時間t以0.01的步長運行了1000次，所以得出來的應該是在前10s內角度的變化，圖像如下：

在前10s內角度的變化

看上去，該函數是一個三角函數，所以當擺角很小的時候，單擺運動可以近似的看作簡諧運動.那麼當擺角逐漸變大時，會發生什麼情況？

釋放角度為3pi/4(135度)
釋放角度為0.99pi(178.2度)
釋放角度為pi(180度)

當釋放角度剛好為 pi 時程序會出現問題，這是因為 sin(pi)=0 ，這就導致角速度一直不會有變化，故角度就不會有變化.所以在最後的圖象上角度總是一個常量 pi .

從上面幾張圖片可以找到一些規律：

1.隨著釋放角度的增大周期也會跟著增大

2.角度的變化始終是周期性的

3.當釋放角度越來越大時，簡諧運動的近似會越來越不精確.

通過上面的程序我們已經初步知道了如何用數值表示這個函數.下次我會講講如何求這個函數的解析解.

參考

[1] 《費恩曼物理學講義》 第一卷 第9章 牛頓的運動學方程


推薦閱讀：
相关文章