插值法(內插與外插)

Interpolation and Extrapolation

 

數學問題探討

當我們在表現資料時,常常會有需要比實際量測點上的值更細密的情況,或者是有需要在範圍外預測其值。

比方說天氣圖的繪製,不論是氣壓或是雨量,都不可能做到處處都有測量站,又例如我們關心一天之中溫度隨時間的變化,但是實際上記錄氣溫的動作可能只是每小時一次,則我們要作一個連續的圖時,就會用到插值法。

插值法的中心議題是:在我們己具備一組表列數(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)

所有的插值法堶掖斨眾瑼熔鶾L於線性插值法,任兩個相鄰的表列點之間必可以拉一條直線把它們連起來,如此在之間的 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 的語法堙A 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 )

 

 

副程式的使用

(1) 注意 polinthunt 或與 locate 的配合。

(2) 使用 polint 時,注意該副程式內己自設上限 NMAX = 10 是最大可用的 n

SUBROUTINE polint(xa,ya,n,x,y,dy)
INTEGER n,NMAX
REAL dy,x,y,xa(n),ya(n)
PARAMETER (NMAX=10) Largest anticipated value of n.

Given arrays xa and ya, each of length n, and given a value x, this routine returns a
value y, and an error estimate dy. If P(x) is the polynomial of degree N ? 1 such that
P(xai) = yai, i = 1, . . . , n, then the returned value y = P(x).

INTEGER i,m,ns
REAL den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)

(3) 要具備使用插值法並呼叫 pgplot 副程式作圖的能力。

 

 

範例副程式(九十七學年度後追加)

以下提供如何使用本節之副程式的範例主程式:

首先第一個例子,這個範例提供 polint 插值的兩種可能用法於同一個主程式中,一個是使用者給一個 x,程式就給一個 y;另一個方法是問範圍及點數,然後程式在此範圍產生等間距的細密插值點(適合作圖用)。

program polint_main

c In subroutine polint, array declaration is used for xa(n) and ya(n),
c for this we need to reclare larger or equal size xa(n) and ya(n) in
c this main program. As for NMAX and c(NAMX) and d(NAMX) in that routine,
c they are taken care locally in that routine alreay with specific size
c (parameter) given when declaring array, therefore we don't need to deal
c with that.

implicit none
integer max_tabu, ngrids_max, n_tabu, m, i, j, k, choice
integer k1, k2, more_points, ngrids
parameter (max_tabu = 100, ngrids_max = 2000)
real xa(max_tabu), ya(max_tabu), x, y, dy, delta_x
real x_plot(ngrids_max), y_plot(ngrids_max),dy_plot(ngrids_max)
real x_min, x_max
character*40 filename1

c Ask filename and the number of tabulated points in the file.

write(*,*) 'Type in the file name of tabulated values (x,y)'
read (*,*) filename1
write(*,*) 'How many points are there in the tabulated set ?'
read (*,*) n_tabu

   if (n_tabu.gt.max_tabu) stop 'Too many tabulated values,
  & please recompile the main program with a larger max_tabu.'

c Open the file and read in the tabulated values, and write to screen
c for double check.

open (unit = 20, file = filename1)

do i = 1, n_tabu
read (20,*) xa(i), ya(i)
write(*,*) xa(i), ya(i)
end do

close (20)

c Ask the order of polynomial to be usedfor interpolation, not higher
c than 10.


1000 continue

write(*,*) 'How many points do you want to form polynomial ?'
write(*,*) 'For eample, 2 points mean linear interpolation.'
read (*,*) m

if (m.gt.10) then
write(*,*) 'Waring : Neumarial Recipe default max is 10.'
write(*,*) 'Please decress your value or set a bigger NMAX'
write(*,*) 'in the subroutine polint.'
write(*,*) 'Your choice (0 : exit) or (1: input again).'
read (*,*) choice
if (choice.eq.0) stop 'Terminated by user.'
if (choice.eq.1) goto 1000
goto 1000
endif

c One can choose to type in x one by one or to generte a continuous set
c for graph plotting.

write(*,*) 'What do you want to do ?'
write(*,*)'0:input x, find y ; 1:generate a set of points (x,y)'
read (*,*) choice

c This part get x one by one, using "locate" to find the right interval.

if (choice.eq.0) then

1002 write(*,*) 'Type in your x:'

read(*,*) x
call locate (xa, n_tabu, x, j)
k = min( max(j-(m-1)/2,1) , n_tabu+1-m )
call polint(xa(k), ya(k), m, x, y, dy)
write(*,*) 'Interpolated value is ', y, ', error is ', dy
write(*,*) 'One more point (yes=1;no=0) ?'
read (*,*) more_points
if (more_points.eq.1) goto 1002
endif

c This part generate a set of x, using "hunt" to find the right interval.
c The resulting (x,y) is written out into a file named by user.

if (choice.eq.1) then

write(*,*) 'Type your range x_initial and x-final:'
read (*,*) x_min, x_max
write(*,*) 'How many points do you want to plot ?'
read (*,*) ngrids
if(ngrids.le.1) stop 'Need at least two grids.'

delta_x = ( x_max - x_min ) / (ngrids - 1)
j = 1
do i=1, ngrids
x = x_min + (i-1)*delta_x
call hunt(xa, n_tabu, x, j)
k = min( max(j-(m-1)/2,1) , n_tabu+1-m )
x_plot(i) = x
write(*,*) 'x is ',x_plot(i),' j is ',j,' k is',k
call polint(xa(k),ya(k),m,x_plot(i),y_plot(i),dy_plot(i))
end do

write(*,*) 'Plrease give a name to your output file:'
read (*,*) filename1
open (unit=30, file=filename1)
do i=1,ngrids
write(30,*) x_plot(i), y_plot(i)
end do
close(30)

write(*,*) 'Info: Error estimate will be written into file
&error.dat'
open (unit=40, file='error.dat')
do i=1,ngrids
write(40,*) x_plot(i), dy_plot(i)
end do
close(40)

end if

end

 

另外,當你要利用插值法來配合繪圖,以下是一個範例程式:

program polint_pgplot

integer ntab_max, nplot_max
parameter (ntab_max = 20, nplot_max = 200)
real x_tab(ntab_max), y_tab(ntab_max), x_min, x_max, delta_x
real x_plot(nplot_max), y_plot(nplot_max), dy_plot(nplot_max)
integer ntab, i, nplot, m, j, pgopen
character*40 filename1

write(*,*)'Please tell me the number of tabulated values:'
read(*,*) ntab

write(*,*) 'Which filename contains the tabulated value ?'
read(*,*) filename1

open(unit=10, file=filename1)
do i=1,ntab
read(10,*) x_tab(i), y_tab(i)
end do
close(10)

write(*,*) 'Type in min and max of x:'
read(*,*) x_min, x_max

write(*,*) 'How many points to draw ?'
read(*,*) nplot

write(*,*) 'How many point to form polynomial ?'
read(*,*) m

delta_x = (x_max - x_min) / (nplot - 1)

do i=1, nplot
x_plot(i) = x_min + (i-1) * delta_x
call hunt(x_tab,ntab,x_plot(i),j)
k = min( max(j-(m-1)/2,1) , ntab+1-m )
call polint(x_tab(k), y_tab(k), m, x_plot(i), y_plot(i),
& dy_plot(i))
end do

do i=1,nplot,10
write(*,*) 'x_plot and y_plot,', x_plot(i), y_plot(i)
end do

if ( pgopen('?') .le. 0 ) stop
call pgenv(2.0, 6.0, 0.0, 40.0, 0, 1)
call pgline(nplot, x_plot, y_plot)
call pgpt(ntab, x_tab, y_tab, 9)
call pgclos

end