北京科技大学数理学院卫宏儒weihr@ustb.edu.cn计算方法第7章插值法插值法是函数逼近的重要方法之一,有着广泛的应用。在生产和实验中,函数f(x)或者其表达式不便于计算复杂或者无表达式而只有函数在给定点的函数值(或其导数值),此时我们希望建立一个简单的而便于计算的函数(x),或为各种离散数据建立连续模型,使其近似的代替f(x),具体有很多种插值法,其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点,常用的插值还有Hermit插值,分段插值和样条插值。求近似函数的方法:由实验或测量的方法得到所求函数y=f(x)在互异点x0,x1,...,xn处的值y0,y1,…,yn,构造一个简单函数p(x)作为函数y=f(x)的近似表达式y=f(x)p(x)使p(x0)=y0,p(x1)=y1,,p(xn)=yn,(a)这类问题称为插值问题。f(x)称为被插值函数,p(x)称为插值函数,x0,x1,...,xn称为插值节点。(a)式称为插值条件。常用的插值函数是多项式。基本概念插值的任务就是由已知的观测点,为物理量(未知量)建立一个简单的、连续的解析模型,以便能根据该模型推测该物理量在非观测点处的特性。估计f(x)在区间[a,b]中某点的值时,当属于包含结点的最小闭区间时,相应的插值称为内插,否则称为外插。在某一逼近函数类中选取的一组线性无关的函数,此时对应的插值函数为:由插值条件确定函数组称为插值基函数。xxx012,,,,nxxxx,0,1,2,,iinpx0011()nnpxcxcxcx(0,1,2,,)icin,0,1,2,,ixin最简单的插值函数是代数多项式Pn(x)=a0+a1x+…+anxn,…...(1)这时插值问题变为:求n次多项式Pn(x),使满足插值条件pn(xi)=yi,,i=0,1,2,…,n,……(2)只要求出Pn(x)的系数a0,a1,…,an即可,为此由插值条件(2)知Pn(x)的系数满足下列n+1个代数方程构成的线性方程组a0+a1x0+…+anx0n=y0a0+a1x1+…+anx1n=y1…………………….a0+a1xn+…+anxnn=yn……(3)而ai(i=0,1,2,…,n)的系数行列式是Vandermonde行列式=……(4)由于xi互异,所以(4)右端不为零,从而方程组(3)的解a0,a1,…an存在且唯一。解出ai(i=0,1,2,…n),Pn(x)就可构造出来了。但遗憾的是方程组(3)是病态方程组,当阶数n越高时,病态越重。为此我们从另一途径来寻求获得Pn(x)的方法----Lagrange插值和Newton插值。xxxxxxxxxxxxnn2nnn1211n0200n10...1..................1...1),...,,V(niijjixx110)(Lagrange插值一、Lagrange插值多项式先从最简单的线性插值(n=1)开始。这时插值问题(2)就是求一次多项式L1(x)=a0+a1x使它满足条件L1(x0)=y0,L1(x1)=y1,令L1(x)=l0(x)y0+l1(x)y1,由于l0(x0)=1,l0(x1)=0,l1(x0)=0,l1(x1)=1.这样l0(x)含有因子x-x1,令l0(x)=λ(x-x1),再利用l0(x0)=1确定其中的系数,结果得到x-x1l0(x)=------------,x0-x1类似的可得到x-x0l1(x)=------------,x1-x0这样。。。(5)l0(x),l1(x)称为以x0,x1为节点的插值基函数。101001011)(yxxxxyxxxxxL线性插值仅仅用两个节点以上的信息,精确度较差。为了提高精确度,我们进一步考察以下三点的插值问题:作二次多项式L2(x)=a0+a1x+a2x2使其满足条件L2(x0)=y0,L2(x1)=y1,L2(x2)=y2令L2(x)=l0(x)y0+l1(x)y1+l2(x)y2。由l0(x0)=1,l0(x1)=0,l0(x2)=0,l1(x0)=0,l1(x1)=1,l1(x2)=0,l2(x0)=0,l2(x1)=0,l2(x2)=1.这样l0(x)含有x-x1,x-x2两个因子,令l0(x)=λ(x-x1)(x-x2),利用l0(x0)=1确定其中的系数λ,得(x-x1)(x-x2)l0(x)=------------------,(x0-x1)(x0-x2)类似的可以得出l1(x),l2(x):(x-x0)(x-x2)(x-x0)(x-x1)l1(x)=-----------------,l2(x)=------------------.(x1-x0)(x1-x2)(x2-x0)(x2-x1)于是(x-x1)(x-x2)(x-x0)(x-x2)(x-x0)(x-x1)L2(x)=-----------------y0+-----------------y1+------------------y2...(6)(x0-x1)(x0-x2)(x1-x0)(x1-x2)(x2-x0)(x2-x1)l0(x),l1(x),l2(x)称为以x0,x1,x2为节点的插值基函数。仿照线性插值和二次插值的办法,进一步讨论一般形式的n次多项式Ln(x)=a0+a1x+a2x2+…+anxn,使其满足Pn(x0)=y0,Pn(x1)=y1,......,Pn(xn)=yn…(7)我们仍从构造插值基函数着手,先对某个固定的下标i,作n次多项式li(x),使其满足条件…(8)容易求得(x-x0)(x-x1)...(x-xi-1)(x-xi+1)...(x-xn)li(x)=-----------------------------------------------=(xi-x0)(xi-x1)...(xi-xi-1)(xi-xi+1)...(xi-xn)0,()1,ijijijlx0njjijijxxxx001110011100()()()()()...()()...()()()()...()()...()()innjjjnjjnnjjjjjjjjjnnnijijjiijxxxxxxxxxxlylLxxxxxyLxxxxxxxxxxxyxx将代入中得....(9)公式(9)就是Lagrange插值多项式,li(x)称为以x0,x1,...,xn为节点的Lagrange插值基函数。二、Lagrange插值的截断误差定理:设Ln(x)是过点x0,x1,x2,…xn的n次插值多项式,,f(n+1)(x)在[a,b]上存在,其中[a,b]是包含点x0,x1,x2,…,xn的任一区间,则对任意给定的x[a,b],总存在一点(a,b)(依赖于x)使…(10)其中,f(n+1)()是f(x)的n+1阶微商在的值。],[)(baCxfn(1)1()()()()()(1)!nnnnRxfxLxxnf101()()()...()nnxxxxxxx证明:记Rn(x)=f(x)-Ln(x)显然Rn(xi)=0,i=0,1,…,n,故可设Rn(x)=K(x)n+1(x)现在[a,b]上任意固定一点x,引进辅助函数g(t)=f(t)-Ln(t)-K(x)n+1(t),(*)则g(t)在[a,b]上具有n阶连续导数,在(a,b)内存在n+1阶导数,在t=x0,x1,…,xn,x诸点处皆等于零,即g(t)在[a,b]中有n+2个零点,由Rolle定理知g'(t)在[a,b]中有n+1个零点,如此反复,最后可推知g(n+1)(t)在[a,b]中有1个零点,,即有g(n+1)()=0,ab.因为n+1(t)是n+1次多项式,n+1(n+1)(t)=(n+1)!,又因为Ln(t)是次数为n的多项式,因此Ln(n+1)(t)=0。这样,由(*)式便有由此得K(x)=f(n+1)()/(n+1)!.代入Rn(x)=K(x)n+1(x),定理得证.(1)(1)(1)(1)()()()()()0nnnnngfLKx上式称为带余项的Lagrange插值公式,只要f(x)具有n+1阶导数,就有上式成立,其余项为特别,当n=1时,取x0=a,x1=b,则有令x1-x0=b-a=h,x=x0+th,0t1则易证,当0t1时,|t(1-t)|的最大值为1/4,1(1)()(1)!()()()nnfnfxpxxabn)()(1)!1()()1(xxRnnfnn)(!2)()(21xfxR22)1()(httx|)(|8|)(|)2(21fhxR应当指出,余项表达式只有在f(x)的高阶导数存在时才能应用。在(a,b)内的具体位置通常不可能给出,如果我们可以求出那么插值多项式pn(x)逼近f(x)的截断误差是…(11)性质:假设x0,x1,…,xn是n+1个互异节点,函数f(x)在这组节点的值f(xk)(k=0,1,…,n)是给定的,那么存在唯一的n次次多项式pn(x)满足pn(xk)=f(xk),k=0,1,…,nMfnnbxax1)1)((max)()!1(11)(xnnMxRnn三、例题:已给sin0.32=0.314567,sin0.34=0.333487,sin0.36=0.352274,用线性插值及抛物插值计算sin0.3367的值并估计截断误差。解:由题意取x0=0.32,y0=0.314567,x1=0.34,y1=0.333487,x2=0.36,y2=0.352274。用线性插值及抛物插值计算,取x0=0.32及x1=0.34,又由公式得y1-y0sin0.3367L1(0.3367)=y0+————(0.3367-x0)x1-x00.01892=0.314567+———(0.0167)=0.330365.0.02其截断误差得其中,因f(x)=sinx,f//(x)=-sinx,可取,于是R1(0.3367)=sin0.3367–L1(0.3367)1/2(0.3335)(0.0167)(0.0033)0.9210–5,若取x1=0.34,x2=0.36为节点,则线性插值为,))((2)(1021xxMRxxx)(//10max2xfxxxM330387.0)0033.0(02.0018787.0333487.0)3367.0(3367.01~0.3367sin112121)(xxxyyyL3335.0)(110)(max2xSinxxxxSinM其截断误差为,其中于是用抛物插值计算sin0.3367时,可得))((2)(~2121xxMxxxR3523.0)(//102maxxfMxxx51036.1)0233.0)(0023.0)(3523.0(21)3367.0(~3367.0sin)3367.0(~11LR330374.00008.0105511.0352274.00004.01089.3333487.00008.0107689.0314567.0)3367.0())(())(())(())(())(())((3367.0sin4442120210221012012010210yLxxxxxxyxxxxxxyxxxxxxxxxxxx这个结果与六位有效数字的正弦函数表完全一样,这说明查表时用二次插值精度已相当高了。其截断误差得其中于是828.0)(0cos///20max3xxMfxxx62210178.0)0233.0)(033.0)(0167.0)(