第三週上課綱要
在本週三小時半的課程中，我們將陸續完成下列幾個單元，每單元約半小時至一小時。
一、Fortran 程式寫作與編譯複習（30 分鐘）
複習寫可計算的程式
如：從 1 加到 N 的程式
複習如何用編譯器做出可執行檔
如：gfortran my_prog.f 或 gfortran -o my_prog.x my_prog.f
二、PGPLOT 繪圖副程式庫的詳細介紹
PGPLOT 簡介（一小時）
示範 pg_test1.f
( gfortran -o pg_test1.x pg_test1.f -L /usr/X11R6/lib -lX11 -L /usr/local/pgplot -lpgplot )
補充資料：Xmgr 互動式繪圖軟體 （30 分鐘）
三、PGPLOT 繪圖之實習（一小時至一小時半）
用打點的方式，畫出三角函數 y = sin(x) 的圖形（提示：call PGPT ）
用連線的方式，畫出 y = x2 sin(x) 的圖形（提示：call PGLINE ）
（第三週實際只做到這裏）
四、 以 PGPLOT 輸出動畫的技巧介紹與實習
換顏色的指令
範例程式閱讀 planet.f
範例程式閱讀 md_pgplot.f
五、電腦模擬實驗：咖啡冷卻問題
問題：先加奶精或是後加奶精咖啡冷得快？
物體量：溫度隨時間的變化
物理定律（數學模型）：牛頓冷卻定律（溫度變化率與溫差呈正比）
dT / dt = - r (T - Ts)
其中溫度 T 是時間 t 的函數，即 T(t) ，Ts 是外在環境的溫度（室溫），
冷卻係數 r 如何定出？
以下為實測值，當時室溫 17 度 C：
時刻（分鐘） 黑咖啡溫度(度C) 白咖啡溫度(度C) 0 82.3 68.8 2 78.5 64.8 4 74.3 62.1 6 70.7 59.9 8 67.6 57.7 10 65.0 55.9 12 62.5 53.9 14 60.1 52.3 16 58.1 50.8 18 56.1 49.5 20 53.3 48.1 22 52.8 46.8 24 51.2 45.9 26 49.9 44.8 28 48.6 43.7 30 47.2 42.6 32 46.1 41.7 34 45.0 40.8 36 43.9 39.9 38 43.0 39.3 40 41.9 38.6 42 41.0 37.7 44 40.1 37.0 46 39.5 36.4 白咖啡與黑咖啡是否都滿足牛頓冷卻定律所描述的呢？如何確認？如何做下一步？
演算法：
Euler 演算法
對於一個簡單的微分方程 dy / dx = f(x,y)，在 Dx 很小時可當作 Dy / Dx = f(x,y)，故 Dy = f(x,y) Dx
又當 Dx 是等間隔 xn - xn-1 = ... = x3 -x2 = x2 -x1 = Dx 時
Dy = yn+1 - yn = f(xn,yn) Dx
這裏的
yn+1 = yn - f(xn,yn) Dx
就是有名的 Euler（奧依勒）演算法
程式流程：給定冷卻係數、咖啡初始溫度、室溫，畫出 T(t) 的圖
(0) 開始
(1) 讀入冷卻係數、咖啡初始溫度、室溫、小間隔時間設定、希望模擬之時間
(2) 在迴圈中，以 Euler 法預測下一小段時間後的溫度，並畫圖，直到迴圈跑完
(3) 結束
範例程式參考
分析與討論
我們怎麼知道程式寫得對或錯？
取 Dt 大一點或小一點有什麼不同？
假設加奶精總是把溫熱咖啡的溫度再降五度，先加奶精過是後加會比較快喝到咖啡？
變化率與當時之量有關的現象，還有什麼例子？
六、電腦模擬實驗：自由落體問題
問題：想知道炮彈打到那裏
物理量：軌跡隨時間的變化 x(t), y(t)
物理定律（數學模型）：f = m a ，在此 f 純粹來自於重力 （粗體字代表是向量）
粒子在地表附近所受的重力是 g （g = GmM/R），g = 9.8 牛頓/米
f = ma = m ( ax ex + ay ey ) ，其中 ax = dvx/dt 、 ay = dvy/dt
只有重力時，f = fy ey = g m ey ，也就是說 fy = g m
演算法：
利用 Euler 演算法：
位置 xn+1 = xn + vx,n Dt ； yn+1 = yn + vy,n Dt
速度 vx,n+1 = vx,n + ax,n Dt ； vy,n+1 = vy,n + ay,n Dt
進階補充：對計算穩定性及精確度，有更好的選擇
Euler-Richarson 演算法， 先求出
ax,n = fx(rn,vn,t)/m
vx,mid = vx,n + (1/2) ax,n Dt
xmid = xn + (1/2) vx,n Dt
ax,mid = fx( xmid, vx,mid, t+(1/2)Dt ) / m
然後才做
vx,n+1 = vx,n + ax,mid Dt
xn+1 = xn + vx,mid Dt
註：為什麼 Euler-Richarson 演算法 比 Euler 的好，可參考 Gould 課本上的推導： p1、p2
程式流程：
(0) 開始，宣告變數
(1) 問質量、初速度（含初速率與角度）
(2) 建立位置 x 與 y、速度 vx 與 vy 的初始值
(3) 用演算法求下一個小時段 Dt 後的位置與速度
(4) 畫出粒子位置並檢查高度是否己達地面，未達地面則回到 (3)
(5) 問要不要再射一次，若要則回到 (1)
(6) 結束
做不出來嗎？參考範例程式
動動腦、動動手，自己再改造（或回家做）：
1. 畫個大砲把炮管的角度表現出來
2. 以火藥包裝填數表現炮彈動能，轉換成初速率
3. 物體被炸破片的彈出（配合亂數產生器 ran3.f）
分析與討論：
想預測炮彈打到那裏，主要是基於那一條物理定律？
什麼樣的先決（初始）條件決定了炮彈的軌跡（和落點）？
同樣的火藥包充填數，怎麼打炮彈才飛得最遠？
七、電腦模擬實驗：受空氣阻力的落體
問題：在有空氣阻力的情況下，炮彈軌跡變得如何？
物理量：一樣是軌跡隨時間的變化 x(t), y(t)
物理定律（數學模型）：f = m a，受力 f 同時包含了重力及空氣阻力
空氣阻力非常複雜（因此飛機或高速火車的阻力設計需要靠風洞實驗或電腦模擬來進行），基本上它是與速率有關，速率越高阻力越大，靜止的物體則沒有空氣阻力。可見它是一個與速率有關的量，常見的簡化公式，有隨速率一次方呈正比的，也有隨兩次方呈正比的，方向是速度反向（負號）的
Fd = - k1 v
Fd' = - k2 v2
由於有 x 與 y 兩個分量，假設速度與水平線的夾角為 i
Fd,x = Fd cos i
Fd,y = Fd sin i
因此，在有關力之部分的公式，
fx = Fd cos i
fy = - m g + Fd sin i
其他的部分，與上一個議題相同
新議題：如何求得室氣阻力係數 k？
利用終端速度，即垂直落下過程中重力與阻力達平衡，使合力為零不再加速，而呈等速運動時的速度。
合力 = 0 = - m g + Fd = - m g + k1 vy2
即 k1 = vy2 / (m g)，此終端速度可藉由實驗量得
一個粉筆頭大小的石子其終端速度大約是每秒三十公尺，另外，一個重 0.254 克 、半徑 0.254 公分的保利龍球有以下的實測數據（引用 Physics Teacher 24, 153 (1986) )
時刻（秒） 位置（公尺） -0.132 0.0 0.0 0.075 0.1 0.260 0.2 0.525 0.3 0.870 0.4 1.27 0.5 1.73 0.6 2.23 0.7 2.77 0.8 3.35 想想看，我們要怎麼從上列的數據求得空氣阻力係數 k？
演算法：
同樣採用 Euler 或 Euler-Richardson 演算法
程式設計：
與上一個議題相同
參考範例程式
分析與討論：
以同樣的力道向上拋一個小石子，在有空氣阻力與無室空氣阻力的兩種情況比較下，那一種會在空中停留較久？
在有空氣阻力的情況下，拋射物的角度要如何設定其距離才會最遠？
進階閱讀：
Gould and Tobochnik, An Introduction to Computer Simulation Methods -- Applications to Physical Systems, Addison Wesley (1996) Chapter 2, Chapter 3