大角度單擺的運動近似(3):運動的數值模擬 第一篇:Bazinga:大角度單擺的運動近似(1):小角度近似第二篇:Bazinga:大角度單擺的運動近似(2):周期的解與橢圓積分 上次我們已經成功解決了周期的表達式以及周期近似,接下來我們把注意力集中到微分方程 的解 上來.此時我們不使用小角度近似,並試圖求該微分方程的數值解.在往後的幾篇文章中我會把自己的級數解法更新上來.在本文會先用一種我臨時受到啟發的解法來解決該問題.靈感出自《費曼物理講義》. 首先對該物理問題所需要的數學基礎做一個鋪墊. 圖片來自www.bing.com本文的核心數學思想是線性近似.對於一個已知函數 ,由自變數在 處增加引起的增量 ,此增量可以表示為 在 時高階無窮小量 ,這時可以定義微分了.不過這裡我們採用近似的記法,即當 很小時, 這個近似模型可以很好的運用在運動學中,比如在本文的情況下 , 所以可以近似的說,當我們的步長 足夠小的時候, 注意在這裡我們定義 ,這說明了對於任意時刻的角度位移,只要知道它上一刻的角度位移和此刻的角速度就可以得到;而此刻的角速度又由上一刻的角速度和上一刻的角加速度給出.此時牛頓第二定律的微分方程給了我們便利,因為它成功的把角加速度和角度位移聯繫了起來: 這也就是說 故只要我們知道了前一時刻的角加速度和角度位移,就可以表示出這一時刻的角度位移.此時我們為了化簡計算,令常數 .同時為了獲得較高的精確度,我們讓 . 我們還需要給出初始條件.考慮到我們研究的是單擺在大角度情況下的運動,那麼我們令 時 ,而角速度自然是 .為了畫出角度位移隨時間變化的圖像,我讓我同學 @天琴心弦 (編程大佬)幫我寫了一段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度)當釋放角度剛好為 時程序會出現問題,這是因為 ,這就導致角速度一直不會有變化,故角度就不會有變化.所以在最後的圖象上角度總是一個常量 .從上面幾張圖片可以找到一些規律:1.隨著釋放角度的增大周期也會跟著增大2.角度的變化始終是周期性的3.當釋放角度越來越大時,簡諧運動的近似會越來越不精確.通過上面的程序我們已經初步知道了如何用數值表示這個函數.下次我會講講如何求這個函數的解析解. 參考[1] 《費恩曼物理學講義》 第一卷 第9章 牛頓的運動學方程 推薦閱讀: 相关文章 {{#data}} {{title}} {{/data}}