插值法（內插與外插）
發表： 2012 年 02 月 17 日 Filed under: 基本數理 發表留言http://163.13.111.54/numerical_methods/nm_units/interpolation_n_extrapolation_intro_n_polynomial.htm
數學問題探討
當我們在表現資料時，常常會有需要比實際量測點上的值更細密的情況，或者是有需要在範圍外預測其值。
比方說天氣圖的繪製，不論是氣壓或是雨量，都不可能做到處處都有測量站，又例如我們關心一天之中溫度隨時間的變化，但是實際上記錄氣溫的動作可能只是每小時一次，則我們要作一個連續的圖時，就會用到插值法。
插值法的中心議題是：在我們己具備一組表列數（tabulated value）的情況下， 如何得出沒被定到之區域的值。
什麼樣的函數才能被插值，這是數學上討論的間問題，參見課文 p.99 第四段。然而，我們會要用到插值法的場合往往都不知道描述對象背後的函數是什麼形式（但相信其有連續的本質），因此我們也只能盡力求真實。
使用插值法所建立的函數，在表列點上一定要重現原本給定的表列值，否則就不是插值法而是函數近似或曲線擬合的間問題了，它們是不一樣的。
插值的作法，很直觀地來講，就是，(1) 先從表列值來獲得函數 f(x)，再 (2) 用函數 f(x) 求出我們所要的任何特定 x 之 f(x) 函數值。然而，比較精密且系統化的數值方法卻不是用這兩個步驟來進行插值，原因是前述兩階段方法對於插值的精密度並沒有控制，效率較差，也比較會有進位誤差。一般在做插值法，是從欲插值點 x 附近的幾個表列點 xi 開始，建立插值函數 f(x)，並且也試著網羅更多表列點來插值，看隨著項數變多誤差會不會變小，如此找出最適合的函數 f(x)。
我們會比較希望演算法在從表列值建立插值用函數時，也能提供誤差分析以供我們或程式來判斷。畢竟可用的插值函數 f(x) 並非唯一，而即便是己設定了採用一種方法，如多項式法，也會有該使用多少項才最恰當的問題。
建立插值函數所需之鄰近表列值個數，我們稱之為插值法的 order（階），較高階未必保證得到較合理的插值，這點在多項式插值法尤其如此，要小心注意。詳見課文中之例圖
上兩圖實線都是原現象背後的真正值，短虛線代表低階多項式插值結果，長虛線代表高階多項式插值結果。明顯可見，case (a) 高階者較準確，而 case (b) 則是低階較準確。
線性插值法（Linear Interpolation）
所有的插值法裏面最簡單的莫過於線性插值法，任兩個相鄰的表列點之間必可以拉一條直線把它們連起來，如此在之間的 x 值就都有線性函數 y(x) 可以對應到，利用直線上的斜率必為固定值的特性，其公式是（以 (x1,y1)、(x2,y2) 為兩個相鄰的表列點為例）：
(y – y1) / (x – x1) = (y2 – y1) / (x2 -x1)
經整理後得
y = [ (y2 -y1) / (x2 – x1) ] (x – x1) + y1
注意等號的右邊全是 x 與常數，我們因此有了 y(x) 的明確公式可用。
我要求大家對於線性插值法這種較簡單的插值法，應該要能在不看參考資料的情況下做出，即自行把式子寫下來，並且把程式寫出來。
多項式插值法
大家都知道兩點唯一決定一條直線（不轉彎）、三點唯一決定一條二次曲線（會轉一次彎）、四點唯一決定一條三次曲線（會轉兩次彎，有反曲點），等等。這些曲線都是以多項式的形式（變數出現時，些是整數次方）。
一個 n – 1 次曲線的多項式雖有像 y = a(n-1)x(n-1) + a(n-2)x(n-2) + …. + a1x +a0 這樣的通式可以表示出，但必須代入 n 個表列值才能定出 an-1 至 a0 那 n 個係數，一下子不易看出。
數學上有一個 Lagrange 多項式公式，它可以由 n 對 (x,y) 值唯一決定 n-1 階多項式，且公式非常好記， 如課文中的式 (3.1.1)
記法參考：寫下每一項都有係數 yi，分母全是 (xi – xj)，其中 i ≠ j，分子則全部都是 (x – xj)，一樣是 i ≠ j ，這樣的項共有 N 個相加。我們的式子最高項次一定是 xN-1，符合 (x – xj) 連乘，並且當 x 恰為某一個 tabulated point xi 時，yi 會因分子分母一模一樣抵消成 1 而留下來，而其他具 yj 的項則因其分子必定有 (x – xi) 而成為零，因此 y = yi，符合表列值的定義。
有了上面 Lagrange polynomial 的定義，拿來寫個程式，對任給的 x 求插值，也並非不可以，只是這樣就少了誤差分析及精密度控制的動作。
真正有用的演算法，是 Neville’s 演算法。它定出 Pi 為原本表列值 yi，而 Pi(i+1) 則代表一般表列點是 xi 與 xi+1 所構成的一階 Lagrange 多項式。
Neville’s algorithm 厲害的地方，是發現由低階多項式，可以經由組合而系統化地得出更高階多項式（例如自 P12 及 P23 得出 P123），根據下列關係式：
另外，不同級數間的大小差異，也方便地定義了C 與 D 如何從 P 得到：
進一步推導，C 與 D 更可以從前一級的 C 與 D 得到
如此，我們再回去看上面那個橫向金字塔，比方說想要得 x 的 P1234，可由 低級數的 P 出發，透過 C 或 D 的求得並附加到 P 上，來獲得較高級數的 P（即在 橫向金字塔 中由左而右）越做越精密， 並且可一路上追蹤不同階數多項式之間的誤差度。
課本提供的副程式是 polint，注意你固然可以有多少個表列點就多少個全用（較危險），你也可以只選其中 m 個（比方說 4 ）來用，只選其中 4 個的呼叫副程式方法，以 tabulated values 是 18 組 xx 與 yy 為例，要用到點 15 ~ 18 的作法是：
call polint( xx(15), yy(15), 4, x, y, dy )
而有別於全部 18 組都使用的
call polint( xx, yy, 18, x, y, dy )
提醒大家，在 Fortran 的語法裏， call polint( xx, yy, 18, x, y, dy ) 就是 call polint( xx(1), yy(1), 18, x, y, dy ) 的意思，即呼叫副程式所傳的引數是陣列或變數的起始位置。
如何搜尋有序的表
我們前面已經建立了從兩個點唯一決定一條直線的線性插值法， 那麼在已知一系列表列點的情況下，被要求要插值某 x 點上求 y，自然我們必需取用 xi < x < xi+1 的那兩個點 ( xi , yi ) 及 ( xi+1 , yi+1 ) 來做線性插值。簡單的說，現在的問題是，給定 x，如何找到 i ？
我們可以想像，若把程式寫成從最小或從最大的表列點開始與 x 比對，萬一 x 值離那端很遠就會沒有效率 。課本提供了二分法 (bisection) 的方法，首先先判斷 x 有沒有小於 x1 或 大於 xN，若確定 x 在其中則拿其中間點 xN/2 （若 N 非偶數就用 N+1）或來與 x 比，判斷出 x 是在 x1 與 xN/2 之間或是 xN/2 與 xN 之間，然後重覆策略，每次都是取用新上下限的中間點去搜尋 。課本提供 locate 副程式給讀者作二分法搜尋，詳見之。 以下圖例：
如果我們要做一系列相鄰插值點的插值，比方說要做圖產生圖點用，則每個新點會鄰近於舊點。若每次都由最大範圍的上下限開始用二分法，則會花很多冤枉工。有效的策略應是，每次找表列點時，從上一次獲得的表列點開始，如此有最大的機會命中，若小於 x，則以次兩倍步幅變大跳躍向右尋找，直到找到的點比 x 大，再改用 bisection，這樣較有效率。 課本提供 hunt 副程式做上述工作，詳見之。 以下圖例：
最後，還有一個問題：雖然用 locate 或 hunt 可以由 x 得 j，其中 j 與 j+1 表列點會把 x 包夾住，但像是多項式插值法，若我們一次要用（一共是 N 個點中的）m 個點（m 比方說是 4），能使用 j-1、j、j+1、j+2 當然很好，但若 j 太靠近兩側，例如 j 是 1 或 j+1 是 N 的話，j 就不能選作是那 m 個點的中間點了，這樣要如何處理？ 答案：使用下列指令，針對 j、m、N 去作運算，則會把邊邊調到剛剛好，為什麼？可自己想一想。
k = min( max(j-(m-1)/2,1) , N+1-m )
呼叫副程式時，像這樣
call polint( xx(k), yy(k), m, x, y, dy )
[C語言數值分析] matrix 矩陣運算 – 問題簡述 / library
發表： 2012 年 02 月 17 日 Filed under: 基本數理 發表留言http://edisonx.pixnet.net/blog/post/83614783
開始前…
在自己下手刻 matrix library 前，請細思：
要自己刻的原因是什麼？因工作需求？因作業需求？還是純粹興趣？
思考完後，或許不需要再看此系列文，反之要做的事情可能是，
去評估一套適合解決手邊問題之 library。
不易解之問題
矩陣問題在線性代數裡佔了不少篇幅，也有不少問題以矩陣模式解之較為合適。
然而在計算機領域，實際解決矩陣問題時，有不少實際狀況必須考量，
較為有名的幾個問題，其中在 CS 領域裡筆者認為「不算好解」的會在後面打 *
另在 (7), (8) 是實作面較為麻煩之事。
(1) 乘法效率 *
如何低於 O(n^3)
(2) 乘法效率 *
一串 matrix 相乘 (如 ABCD) ，相乘次序為何 (如 AB(CD) , A(BC)D )，才能使得在結果正確情況下，乘法次數最少、效率最高？
(3) ill-condition 問題
消除 ill-condition ，這算是較簡單解決之問題。
(4) 消去法效率
消去法種類必須如何低於 O(n^3)
(5) 大型矩陣問題 *
一些較大問題變數個數不小，若在記憶體容量不足之情況下，是否有辦法將大型矩陣切成數個小型矩陣再進行求解動作。
(6) 稀疏矩陣問題 *
稀疏矩陣在資料結構之表示法並不難，難的是如何使其執行效率高。
(7) 如何包成 friendly use class ?
實際包成 class 後發現主要幾個問題點：該用 shallow copy / deep copy、繼承架構該如何設計 、包成 class 後實作上會有不少 side effect 影響效率 (特別是隱喻之 Constructor 便會耗去不少時間) 等，在考量 friendly、effective、multi-thread 情況下，這問題會愈變愈複雜。
(8) 例外處理
目前絕大多數之數值分析在探討時任何 matrix 問題確實會做些基本之例外處理，但一些問題之例外處理並不容易完全處理掉。
Free library
然而實際上也有不少免費 matrix library 已提供，如
1. Matrix Template Library 4 (MLT4) * (不少人用過覺得不錯)
2. C++ Matrix Template Class Library
3. Matrix Expression Templates (MET)
4. Sparse Lib++ * (少數專門處理稀疏矩陣)
5. CAM Matrix / Vector / Array Classes
6. Netmat
7. High Mat C++
(聲明，筆者沒用過任何 matrix library，據說 MLT4 還不錯)
其他有些 Graphic Library ( OpenGL 還是 OpenCV ) 可能也有提供，
也可能有些 library 效能更佳沒紀錄到，也可能有些是 non-Free，
此處便不一一列舉。
參數式的意義
發表： 2012 年 02 月 12 日 Filed under: 基本數理 發表留言參數式是由"一個點"加上一個"方向向量"所組成,即形成直線的參數式
參數式就是表示直線,無論在空間或在平面
參數式是表示那個方向上的所有的點,點即構成直線
用參數式表示點,即那條線上所有可能的位置
ex:(1+t,2+t,3+t) 是假設那條線上的某一點
Like 669699 