Lecture 2 Vpython動力學模擬 - F=ma的數值解

建國高中特色選修課程 - 物理現象的程式設計與模擬
作者：曾靖夫
日期：2018/7/5

一、牛頓第二運動定律與迴圈

在牛頓力學中，牛頓第二運動定律是我們處理所有運動都會使用的定律。原因是你只要知道初始條件，當受力決定了速度就決定了（$\Delta \vec{v} = \vec{a}\ \Delta t$），而速度決定了位置就決定了（$\Delta \vec{r} = \vec{v}\ \Delta t$）。因此你只要知道物體受了什麼力，就可以預測未來每一瞬間的運動狀態。你也可以說這是力學問題的暴力解法，因為即使作用力再複雜都一定可以算出正確的物體運動狀態，只是不是用手算而是讓電腦幫我們算，即是在科學計算所說的「數值解」。

這裡簡單介紹一下牛頓第二運動定律的基本內容

在這個基礎下，使用迴圈去模擬每個時間間隔所發生的變化是最適當的，只要設定的時間間隔夠短，模擬出來的結果就近似為每一時刻的瞬時狀態。以下簡單分析每個迴圈產生的變化：

每個迴圈所計算的東西，會累加到下一個迴圈中，產生一連串的連續動態的變化，這就是我們進行動力學模擬的重要基礎。

二、自由落體運動的模擬（均勻重力場）

這裡我們馬上來分析最熟悉的變速度運動－等加速度運動，這裡以均勻重力場中不受阻力的拋體運動為例，物體只受到固定的重力作用，因此加速度也必為定值。首先我們分析初速度為零的自由落體運動，參考程式如下：

Exapmle 1 : 等加速度運動－自由落體

from vpython import *

g = 9.8            #重力加速度 9.8 m/s^2
size = 0.5         #球半徑 0.5 m
height = 15.0      #球初始高度 15 m
m = 1.0            #球質量1kg

scene = canvas(width=600, height=600,x=0, y=0, center = vector(0,height/2,0)) #設定畫面
floor = box(length=20, height=0.01, width=10, color=color.green)  	#畫地板
ball = sphere(radius = size, color=color.yellow, make_trail= True, trail_type="points", interval=100) 	#畫球

ball.pos = vector(0, height, 0)    #球初始位置
ball.v = vector(0, 0, 0)           #球初速 

dt = 0.001	#時間間隔 0.001 秒
t = 0.0		#模擬初始時間為0秒

while ball.pos.y > size:    #模擬直到球落地 即y=球半徑
    rate(1/dt)    #每一秒跑 1000 次
    t = t + dt    #計時器
    
    ball.a = vector(0,-g,0)/m            #球的加速度
    ball.v = ball.v + ball.a*dt          #球的末速度 = 前一刻速度 + 加速度*時間間隔
    ball.pos = ball.pos + ball.v * dt    #球的末位置 = 前一刻位置 + 速度*時間間隔
    
print (t, ball.v)

執行後物件視窗的畫面如下：

接下來如果我們要使小球碰到地面之後反彈，則需要引入程式撰寫中很常用到的條件判斷結構－if和else。這裡我們假設球碰到地面的時候，以相同度速率彈起，如果整個過程沒有受到任何阻力，小球應該會達到與釋放時相同的高度，在程式中我們很容易驗證這個觀念。在程式中只須增加一點點while迴圈的內容，程式範例如下：

while True:	#無窮迴圈
    rate(1000) 

    ball.v.y += -ball.a*dt    #此寫法意義上與ball.v = ball.v + ball.a*dt完全相同
    ball.pos += ball.v*dt     #此寫法意義上與ball.pos = ball.pos + ball.v*dt完全相同

    if ball.pos.y <= size and ball.v.y < 0:    #條件：球心高度小於球半徑且速度沿-y軸
        ball.v.y = - ball.v.y    #條件成立則球的速度加一負號表示反彈

到目前為止，相信同學已經知道如何寫出在均勻重力場中的等加速度運動，我們只要再加入各種初始條件就可以模擬出各種拋體運動。比如說，你只需要將while迴圈上方的初始速度設定為ball.v = vector(5,5,0)，就會發現球將以仰角 $45^{\circ}$、速率 $5\sqrt{2}$ 斜向拋出。你也可以用三角函數來描述拋射方向，將球度初始速度設定為 $\left ( v_0\ {\rm sin}\theta , v_0\ {\rm cos}\theta, 0 \right )$，將發現球以初速 $v_0$、仰角 $\theta$ 沿 $x-y$ 平面斜向拋出，你也將發現拋體運動是一種很典型的平面運動。

以上我們討論了無空氣阻力的拋體運動，但在現實生活中，這樣的運動太過理想。其實我們也知道空氣阻力對物體運動的影響其實不小，這裡將試著引入空氣阻力。一般而言，在物體運動速度不大的情況下，物體受空氣阻力的大小與其速率成正比，阻力方向與速度方向相反，數學式可以表示為：

課堂作業 2-1

請完成拋體運動的模擬，使球碰到地板會發生反彈，並使用指令描繪出球的運動軌跡。接著請再新增一個球並加入空氣阻力的影響，保留沒有空氣阻力拋體的軌跡與之對照，空氣阻力常數 $k$ 可以自行設定。

執行成果如上圖所示，黃色球為無空氣阻力，紅色為加入空氣阻力的結果。

三、拋體運動的運動方程式

在這個段落我們還是簡單列出拋體運動的數學方程式，但並不仔細證明推導，目的在於讓同學理解程式工具的重要性。利用程式我們只要用很簡單的觀念將「力」寫入程式中，而不需要經過繁雜的數學推導，就可將物體的運動狀態完整模擬。

如果拋體運動沒有受到空氣阻力，只受到均勻向下的重力作用，我們可以將物體的運動分解為兩個方向，一為沿著水平面的方向，二為沿著鉛直方向。在水平方向上物體並沒有受力，因此在這個方向物體的運動應為等速度運動；而在鉛直方向物體受到固定的重力，所以在這個方向為等加速度運動。

如上圖所示，我們分別由上方和側面照光，就可以透過影子很清楚的呈現它在兩個方向上的運動狀態，按這個概念速度速度向量隨時間的關係為：

進階作業 2-1

請利用程式驗證今天在程式中加入空氣阻力的拋體運動模擬與使用牛頓力學解出的運動方程式完全一致。請在視窗中留下三個球與其運動軌跡，如下圖所示，黃色為無阻力的理想拋體，紅色為加入空氣阻力的拋體，綠色為上方方程式所畫出來的軌跡。這邊我們要求執行流程為，黃球和紅球同時拋出，當紅球和黃球都落地就讓它們直接停在地面上，之後再拋出綠球，將發現綠球會貼著黃色軌跡走，確認「數值解」和「理論解」是吻合的。

四、三角函數簡介

在物理學中很常利用三角函數來表示空間向量在各個方向上的分量。而在數學上三角函數的定義非常單純，只是直角三角形的邊角關係，以下我們對照直角座標系，給出三種三角函數的定義：

1. 餘弦函數：$\ {\rm cos} \theta = \dfrac{x}{r}$ （投影到 $x$ 軸）
2. 正弦函數：$\ {\rm sin} \theta = \dfrac{y}{r}$ （投影到 $y$ 軸）
3. 正切函數：$\ {\rm tan} \theta = \dfrac{y}{x}$ （斜率）

這裡的角度 $\theta$ 一般可以用度（$\rm deg$）來作為單位，但是在Python中使用更有意義的單位稱為弧度（$\rm rad$）。弧度的定義為，圓心角所對應張開的弧長除以半徑，上圖是 $1 \rm rad$ 的定義因此一個圓周角即為 弧度，因此一般熟知的圓周角 $360\ {\rm deg} = 2\pi\ {\rm rad}$。

同學們知道了三角函數的定義之後，就可以了解在平面上如果有一向量 $\vec{r}$，則此向量可以透過其投影在 $x$ 軸和 $y$ 軸的分量來描述，一般使用直角座標表示法寫為$\left ( r\ {\rm sin}\theta , r\ {\rm cos}\theta \right )$，而在程式中我們就是用相同的方法來描述空間向量，只是更多了 $z$ 方向的空間維度。

從以上的座標圖，你也將發現 $\sqrt{x^2+y^2}=r$，這就是圓方程式，從圓心到圓弧上一點的距離是固定的。所以我們可以很方便的用三角函數去描述一些固定大小，但是不同方向的物理量。例如：在拋體運動中，如果固定拋射速率，只改變不同的拋射角度，則我們就會將初速度向量寫為 $\vec{v_{0}}=\left ( v_0\ {\rm sin}\theta , v_0\ {\rm cos}\theta, 0 \right )$，意思是將物體以初速 $v_{0}$，在 $x-y$ 平面上以仰角為 $\theta$ 拋出。或許同學可以透過迴圈控制計算出仰角為何的拋體，會具有最大的水平射程。同學們可以試著練習看看，利用Python算算看是否和你的經驗吻合。

進階作業 2-2

固定初速 $v_{0}$，拋射角度由 $\theta=3\ {\rm deg}$ 開始，每次完成拋射增加 $3\ {\rm deg}$ ：

1. 請在物件視窗，做出多個拋體運動並留下軌跡，同時也使用label指令在畫面上產生一個標籤顯示每次拋射時的角度，並將水平射程R、飛行時間T對上角度theta的關係，以graph指令呈現函數圖，如下所示：

   label指令請參考 vpython7 - label 或 vpython6 - label。

2. 請在執行視窗留下水平射程數據，並驗證 $\theta=45\ {\rm deg}$ 時具有最大水平射程。

五、等速率圓周運動

同學們了解三角函數的定義後，可以馬上利用 $x$ 和 $y$ 方向的分量合成出圓軌跡，我們只要在迴圈中每一次增加一點小小的角度，迴圈持續的跑下去就會變成是動態的圓周運動。

現在我們做一個回轉半徑 $R=1$ 公尺，角速率大小 $\omega = 2\pi \ {\rm rad/s}$的等速率圓周運動。這裡的角速度大小指的是每單位時間掃過的角度，$2\pi \ {\rm rad/s}$代表每秒鐘繞圓周一圈。

Exapmle 2 : 等速率圓周運動

from vpython import *

size = 0.1     #球的大小
theta = 0.0    #初始角度
R = 1.0        #圓周運動半徑
omega = 2*pi   #角速度大小=單位時間繞過的角度
t = 0.0        #初始時間

scene = canvas(width=500, height=500, center=vector(0,0,0), background=vector(148.0/225,228.0/225,204.0/225))
ball = sphere(radius=size, color=color.blue, make_trail=True, interval = 1, retain = 900)

ball.pos = vector(R,0,0)    #球的初始位置
t = 0.0       #初始時間
dt = 0.001    #時間間隔

while True:
    rate(1/dt)
    theta += omega*dt    #t時刻的角度
    ball.pos = vector(R*cos(theta), R*sin(theta), 0)
    t += dt              #計時器

到這邊我們已經能做出圓周運動的模擬了，再進一步可以利用這個圓周運動來設計計時功能，以下提供一個簡單的視覺化計時器範例。

Example 3 : 三點記錄法 - 讓程式自動錄下依順序三個相鄰時間點

pre_theta = theta

while True:
    rate(1/dt)
    t = t + dt
    pre_pre_theta = pre_theta    #前前時刻的角度
    pre_theta = theta            #前一時刻的角度
    theta = theta + omega*dt     #現在時刻的角度

    if pre_theta%(2*pi) > pre_pre_theta%(2*pi) and pre_theta%(2*pi) > theta%(2*pi):
        print ('period = ', t)
        #print pre_pre_theta%(2*pi), pre_theta%(2*pi), theta%(2*pi)
        t = 0    #這行不寫t就不會歸零，而會印出t = 1, 2, 3, 4, ...

請同學「有腦的」複製 Example 3 的 code 貼到 Example 2 中，使物體回到原點的時候，都在執行視窗印出記錄數值。

三點記錄法簡化版

同學可以將下面這小段程式執行看看，就能清楚電腦是怎麼錄下依序的三個時刻。

import time  #python內建時鐘

t = 0        #初始時刻
dt = 1       #時間間隔
pre_t = t    #只是為了防出錯，可以隨便設，第2個迴圈起就正常了

while True:
    time.sleep(1)        #每隔1秒跑一次迴圈
    pre_pre_t = pre_t    #前前時刻
    pre_t = t            #前一時刻
    t = t + dt           #現在時刻

    print ('pre_pre_t =', pre_pre_t, ', pre_t =', pre_t, ', t =', t)

執行結果如下：

課堂作業 2-2：

請完成視覺化的圓周圓周運動計時器，並在執行視窗紀錄每繞一圈所花的時間。執行畫面如下圖所示，物件視窗中當物體跨過原點時，執行視窗就會打印出時間。

進階作業 2-3

延續課堂作業 2-2，請將此圓周運動用 cylinder 指令每隔一小段時間就畫出細線，將此圓周分割成 N 等分。執行畫面如下：

上左圖為圓周運動繞第一圈時，這一圈用來計算週期不劃線，等到回到出發點瞬間，電腦將記錄下繞一圈的週期 period_t，然後再將此時間除以 N，於第二圈開始每隔 period_t/N 用 cylinder 指令畫一條線，如上右圖所示為切成 20 等分的樣子。

你可能會疑惑，週期不是一開始就知道，為什麼要大費周章繞完一圈讓電腦求出週期？在 Example 2 中我們預先設定了角速率為 $2 \pi$，也就是一秒繞一圈，那週期早就知道了不是嗎？

事實上，有很多週期性運動我們一開始並不知道週期，例如：行星受萬有引力繞太陽運行，我們沒法從初速和受力就心算出繞行週期，因此必須要等它繞回原點利用程式裡的計時器把週期測出來。這個技巧會在未來 Lecture 4 克普勒第二定律用上，事實上行星運動是橢圓軌道，你更難一開始就猜到週期究竟是多少。

以下我們提供建議的寫法：

作業繳交：

平台連結：…

此平台由南港高中-高慧君老師設計建置