插值法内插与外插

插值法内插与外插
插值法(内插与外插)
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)
所有的插值法里面最简单的莫过于线性插值法，任两个相邻的表列点之间必可以拉一条直线把它们连起来，如此在之间的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 xxi+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)
副程式的使用
注意polint、hunt或与locate的配合。

使用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 avalue 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)
要具备使用插值法并呼叫pgplot副程式作图的能力。

范例副程式(九十七学年度后追加)
以下提供如何使用本节之副程式的范例主程式：
program polint_mainc In subroutine polint,array declaration is used for xa(n)and ya(n),
c for this we nee
d to reclar
e larger or equal size xa(n)and
ya(n)in cthis 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 nee
d to deal cwith 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 filename1c 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 alarger max_tabu.'
c Open the file an
d read in th
e tabulated values,and write to screen cfor 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 cthan 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 abigger NMAX'
write(*,*)'in the subroutine polint.'
write(*,*)'Your choice(0：exit)or(1：input again).'
read(*,*)choice if(choice.eq.0)stop'Terminaed by user.'
if(choice.eq.1)goto 1000 goto 1000 endifc One can choose to type in xone by one or to generte acontinuous set cfor graph plotting.
write(*,*)'What do you want to do?'
write(*,*)'0：input x,find y；1：generate aset of points(x,y)'
read(*,*)choicec This part get xone by one,using"locate"to find the right interval.
if(choice.eq.0)then1002 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 endifc This part generate aset of x,using"hunt"to find the right interval.
c The resulting(x,y)is written out into afile name
d 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 aname 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(*,*)'Ple ase 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。