內(nèi)插與外插要點

當我們在表現(xiàn)資料時，常常會有需要比實際量測點上的值更細密的情況，或者是有需要在范圍外預測其值。比方說天氣圖的繪制，不論是氣壓或是雨量，都不可能做到處處都有測量站，又例如我們關(guān)心一天之中溫度隨時間的變化，但是實際上記錄氣溫的動作可能只是每小時一次，則我們要作一個連續(xù)的圖時，就會用到插值法。插值法的中心議題是：在我們己具備一組表列數(shù)（tabulatedvalue）的情況下，如何得出沒被定到之區(qū)域的值。我們會要用到插值法的場合往往都不知道描述對象背後的函數(shù)是什麼形式（但相信其有連續(xù)的本質(zhì)），因此我們也只能盡力求真實。使用插值法所建立的函數(shù)，在表列點上一定要重現(xiàn)原本給定的表列值，否則就不是插值法而是函數(shù)近似或曲線擬合的間問題了，它們是不一樣的。插值的作法，很直觀地來講，就是，(1)先從表列值來獲得函數(shù)f(x)，再(2)用函數(shù)f(x)求出我們所要的任何特定x之f(x)函數(shù)值。然而，比較精密且系統(tǒng)化的數(shù)值方法卻不是用這兩個步驟來進行插值，原因是前述兩階段方法對於插值的精密度并沒有控制，效率較差，也比較會有進位誤差。一般在做插值法，是從欲插值點x附近的幾個表列點xi開始，建立插值函數(shù)f(x)，并且也試著網(wǎng)羅更多表列點來插值，看隨著項數(shù)變多誤差會不會變小，如此找出最適合的函數(shù)f(x)。我們會比較希望演算法在從表列值建立插值用函數(shù)時，也能提供誤差分析以供我們或程式來判斷。畢竟可用的插值函數(shù)f(x)并非唯一，而即便是己設定了采用一種方法，如多項式法，也會有該使用多少項才最恰當?shù)膯栴}。建立插值函數(shù)所需之鄰近表列值個數(shù)，我們稱之為插值法的order（階），較高階未必保證得到較合理的插值，這點在多項式插值法尤其如此，要小心注意。詳見課文中之例圖兩圖實線都是原現(xiàn)象背後的真正值，短虛線代表低階多項式插值結(jié)果，長虛線代表高階多項式插值結(jié)果。明顯可見，case(a)高階者較準確，而case(b)則是低階較準確。線性插值法（LinearInterpolation）所有的插值法里面最簡單的莫過於線性插值法，任兩個相鄰的表列點之間必可以拉一條直線把它們連起來，如此在之間的x值就都有線性函數(shù)y(x)可以對應到，利用直線上的斜率必為固定值的特性，其公式是（以(x1,y1)、(x2,y2)為兩個相鄰的表列點為例）：(y-y1)/(x-x1)=(y2-y1)/(x2-x1)經(jīng)整理後得y=[(y2-y1)/(x2-x1)](x-x1)+y1注意等號的右邊全是x與常數(shù)，我們因此有了y(x)的明確公式可用。提醒大家，在Fortran的語法里，callpolint(xx,yy,18,x,y,dy)就是callpolint(xx(1),yy(1),18,x,y,dy)的意思，即呼叫副程式所傳的引數(shù)是陣列或變數(shù)的起始位置。如何搜尋有序的表我們前面已經(jīng)建立了從兩個點唯一決定一條直線的線性插值法，那麼在已知一系列表列點的情況下，被要求要插值某x點上求y，自然我們必需取用xi<x<xi+1的那兩個點(xi,yi)及(xi+1,yi+1)來做線性插值。簡單的說，現(xiàn)在的問題是，給定x，如何找到i？我們可以想像，若把程式寫成從最小或從最大的表列點開始與x比對，萬一x值離那端很遠就會沒有效率。課本提供了二分法(bisection)的方法，首先先判斷x有沒有小於x1或大於xN，若確定x在其中則拿其中間點xN/2（若N非偶數(shù)就用N+1）或來與x比，判斷出x是在x1與xN/2之間或是xN/2與xN之間，然後重覆策略，每次都是取用新上下限的中間點去搜尋。課本提供locate副程式給讀者作二分法搜尋，詳見之。以下圖例：如果我們要做一系列相鄰插值點的插值，比方說要做圖產(chǎn)生圖點用，則每個新點會鄰近於舊點。若每次都由最大范圍的上下限開始用二分法，則會花很多冤枉工。有效的策略應是，每次找表列點時，從上一次獲得的表列點開始，如此有最大的機會命中，若小於x，則以次兩倍步幅變大跳躍向右尋找，直到找到的點比x大，再改用bisection，這樣較有效率。課本提供hunt副程式做上述工作，詳見之。以下圖例：最後，還有一個問題：雖然用locate或hunt可以由x得j，其中j與j+1表列點會把x包夾住，但像是多項式插值法，若我們一次要用（一共是N個點中的）m個點（m比方說是4），能使用j-1、j、j+1、j+2當然很好，但若j太靠近兩側(cè)，例如j是1或j+1是N的話，j就不能選作是那m個點的中間點了，這樣要如何處理？答案：使用下列指令，針對j、m、N去作運算，則會把邊邊調(diào)到剛剛好，為什麼？可自己想一想。k=min(max(j-(m-1)/2,1),N+1-m)呼叫副程式時，像這樣callpolint(xx(k),yy(k),m,x,y,dy)副程式的使用注意polint、hunt或與locate的配合。使用polint時，注意NMAX=10是最大可用的n值SUBROUTINEpolint(xa,ya,n,x,y,dy)INTEGERn,NMAXREALdy,x,y,xa(n),ya(n)PARAMETER(NMAX=10)Largestanticipatedvalueofn.Givenarraysxaandya,eachoflengthn,andgivenavaluex,thisroutinereturnsavaluey,andanerrorestimatedy.IfP(x)isthepolynomialofdegreeN?