内插与外插要点

当我们在表现资料时，常常会有需要比实际量测点上的值更细密的情况，或者是有需要在范围外预测其值。比方说天气图的绘制，不论是气压或是雨量，都不可能做到处处都有测量站，又例如我们关心一天之中温度随时间的变化，但是实际上记录气温的动作可能只是每小时一次，则我们要作一个连续的图时，就会用到插值法。插值法的中心议题是：在我们己具备一组表列数（tabulatedvalue）的情况下，如何得出没被定到之区域的值。我们会要用到插值法的场合往往都不知道描述对象背後的函数是什麽形式（但相信其有连续的本质），因此我们也只能尽力求真实。使用插值法所建立的函数，在表列点上一定要重现原本给定的表列值，否则就不是插值法而是函数近似或曲线拟合的间问题了，它们是不一样的。插值的作法，很直观地来讲，就是，(1)先从表列值来获得函数f(x)，再(2)用函数f(x)求出我们所要的任何特定x之f(x)函数值。然而，比较精密且系统化的数值方法却不是用这两个步骤来进行插值，原因是前述两阶段方法对於插值的精密度并没有控制，效率较差，也比较会有进位误差。一般在做插值法，是从欲插值点x附近的几个表列点xi开始，建立插值函数f(x)，并且也试着网罗更多表列点来插值，看随着项数变多误差会不会变小，如此找出最适合的函数f(x)。我们会比较希望演算法在从表列值建立插值用函数时，也能提供误差分析以供我们或程式来判断。毕竟可用的插值函数f(x)并非唯一，而即便是己设定了采用一种方法，如多项式法，也会有该使用多少项才最恰当的问题。建立插值函数所需之邻近表列值个数，我们称之为插值法的order（阶），较高阶未必保证得到较合理的插值，这点在多项式插值法尤其如此，要小心注意。详见课文中之例图两图实线都是原现象背後的真正值，短虚线代表低阶多项式插值结果，长虚线代表高阶多项式插值结果。明显可见，case(a)高阶者较准确，而case(b)则是低阶较准确。线性插值法（LinearInterpolation）所有的插值法里面最简单的莫过於线性插值法，任两个相邻的表列点之间必可以拉一条直线把它们连起来，如此在之间的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)的明确公式可用。提醒大家，在Fortran的语法里，callpolint(xx,yy,18,x,y,dy)就是callpolint(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)呼叫副程式时，像这样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?