咖啡冷卻問題裏的大問題：定出 r 值
如何定出 r 值？
像牛頓冷卻定律這樣簡單的公式，我們只需要知道冷卻係數 r 就可以進行模擬，但是，如何從一組實驗量測的溫度對時間變化來找出 r 呢？
直觀的想法與作法：
首先，把溫度對時間的實驗數據畫出，再來，試用不同的 r 值把冷卻曲線畫在一起， 一條接一條用不同的顏色，直到你找到能使模擬曲線看起來最接近的 r 值。
小挑戰：你如何讓程式可以允許你不斷輸入新 r 值繪圖，而畫面上又只保留兩三條線以便看清楚？
為了讓大家省卻輸入的麻煩，在此已經將 時刻-黑咖啡溫度-白咖啡溫度 存在一個文字檔中，大家可用 Fortran 開檔 (OPEN) 讀入 (READ) 兩種指令的組合，把數據直接從外部讀進來。範例程式片段如下：
open (unit=10, file='coffee_expt.dat')
do j = 1, j_max
read (10,*) time_expt(j), T_white(j), T_black(j)
end do
close (10)
以上的程式，是利用開啟了陣列的方式，來把存入自外部讀入的數據，由於 Fortran 77 在一開始宣告陣列時就必須先確定陣列的大小，通常有兩種策略來配合進行宣告，分述如下：
第一種，是宣告一個比較大一點陣列，足夠檔案的資料長度所使用，頂多再測試一下迴圈次數是否會超過陣列大小。如下：
real time_expt(100), T_white(100), T_black(100)
...
j_max = 24 （或 read (*,*) j_max ）
if (j_max .gt. 100) then
write(*,*) 'Error, data length is loger than arry size.'
stop
end if
...
do j=1, j_max
第二種常用的方法，是把宣告了整數變數，進一步以 parameter 指令設定變成常數（一旦設定成常數則在程式中不會改變），如下面方式：
integer j_max
parameter (j_max = 24)
real time_expt(j_max), T_white(j_max), T_black(j_max)
...
do j=1, j_max
老師寫的範例程式：
cooling_r_n_expt.f 、 cooling_r_n_expt.f.txt
量測的數據：
較嚴密的作法：
先想什麼叫做接近？那什麼又叫做兩條曲線接近？
跟遠、近有關的物理量是什麼？ 答案就是距離，所謂的接近，就是兩點間距離短。
至於兩條線的情形如何？先想兩種極端的情況，兩條曲線完全重合一模一樣最接近，形狀位置完全不同，就會差比較遠。線是許多點構成的集合，對於模擬出來的曲線 T(t) 與實驗數據 {(tj, Tj), for all j} 而言，我們也可以定義一種距離 d(r) = F(r)(1/2)，其中 F(r) 是以下函數
F(r) = Σj [T(tj) - Tj]2
由於 x(1/2) 是一個遞增函數，所以能使 F(r) 最小的 r 值就會是同一個能使 d(r) 最小的 r 值，我們只要問什麼樣的 r 能使 F(r) 達到極小即可。
接下來，是要寫一段程式來進行 F(r) 函數值的計算。這就要進行以下的動作：
(1) 把實驗值 time_expt(j) 及 T_expt(j) 準備好，
(2) 在給定的 r 值情況下，跑一輪模擬的 T(t)，也就是有一組遠比 j 點數更多點的 time(i) 及 T(i)，在進行迴圈 i 的過程中，要
(3) 研判下一個 j 是否剛好超過，若是，則以 i-1 與 i 兩點來進行所謂的線性插值。插值的意思，就是好比說我們已知第 87 秒是 45.0 度、第 88 秒是 46.0 度，問在 87.75 秒是幾度這種問題，使用線性差值，就是把兩個己知道數據點之間未知的部分直接拉一條直線來代表，這是一種最簡單的取求插值的方法，因為我們的 i 迴圈每一格是推進很小的 Dt 時間，所以線性差值很合理。
小註解 1：線性差值
在已知兩點 (x1,y1)、(x2,y2) 點之間的任意點 x，其 y 值（由於 (x,y) 是在同一條直線上）必定滿足以下關係式
(y2 - y1) / (x2 - x1) = (y -y1) / (x - x1)
整理後可得
y = S (x - x1) + y1
其中 S 是直線的斜率，也就是 (y2 - y1) / (x2 - x1) 。
小註解 2：自變數重合
如果覺得上面的步驟 (3) 太麻煩，在這裏由於實驗數據的量測時間點是有非常規律間距的特珠的情況，如第 0 分鐘、第 2 分鐘、第 4 分鐘等，有剛好的規律跟間距，因此，時間間隔 Dt 可以安排成剛好會有模擬的時刻與量測的時刻點是重合的，如此就不必再借助上述的步驟 (3) 了。
符合以上功能的程式，會問初始及最終的 r，以及 r 的變化要分幾段的範例程式如下：
（現場實作，大家加油）
最嚴謹的作法：
前述的方法藉由作圖目視來決定最好的 r 值，目測一個拋物線（只要所看的範圍夠小，平滑函數枉極值發生之處一定是拋物線的形狀）底部的位置雖然已經比在一堆線之由找最相接近的要容易得多，但仍不是全自動的方法。最嚴謹的作法，直接用數值演算法的技巧將函數 F 最小化找 r，方法可在數值計算方法的教科書找到，在此暫不作深入介紹。
定出 r 值以外的問題：模型與數據的吻合度探討
前述的作法固然卻是明確地足以給出 r 值，但光得到 r 卻沒有讓我們知道所採用的數學模型（物理定律）本身它對現象的描述是很吻合或不太吻合實驗數據（千萬不要忘了實驗數據本身也含有誤差）。檢視上述的函數 F(r) 的大小，或者是是 F(r) 的開根號 d 值，提供我們一個量度的指標。
事實上，F(r) 在形式上很像是統計學上是變異數的定義，而 d(r) 則像是標準差。在統計的意義上，標準差越小，代表數據偏離與散亂的程度越低，而在我們的例子裏，則是意味著數學模型越符合現象的觀察。
延伸思考：
r 值是否為常數？
在本單元中，冷卻係數 r 的值在冷卻定律公式中是被視作為一個常數。如果它不是一個常數，是否對擬合（重現）實驗數據的一致度上會更提昇？
在此所謂不是常數，可以有兩種可能，一是 r 可隨時刻變化，或是可隨溫度（或溫差）變化。其中前者的可能性將導致時間對稱性的問題，除非狀況特殊（例如液體成份隨時間變化而逐漸變質），否則還是以第二種情況，即 r 本身隨溫度（差）的變化仍有相關，是我們較有興趣探討的可能方向。
如果牛頓冷卻定律中的 r 仍是溫度（差）的函數，我們可以怎樣來處理這個數學模型的來解呢？
基於原有之 r 為溫度（差）常數的模型已經能大致地重現出不錯的冷卻曲線，我們可以假設常數加上一個相當微小的溫度一次項可能把數據重現
如此，則數學模型就變成為含有溫度（差）的一次方及平方項，我們可以再次使用最小平方法作曲線擬合，使計算值
除了最小化差的平方和之外，用邊界條件法（如射擊法）也可以得到最佳的兩個待定參數