1、2. 插值法 在生产和实验中，常常需要根据一张表格表示的函在生产和实验中，常常需要根据一张表格表示的函数推算该表中没有的函数值数推算该表中没有的函数值.解决此类问题的简单途径之解决此类问题的简单途径之一利用插值法。一利用插值法。 插值在数学发展史上是一个老问题，它是和插值在数学发展史上是一个老问题，它是和Gauss, Lagrange, Newton等在著名数学家连在一起的。它最初等在著名数学家连在一起的。它最初来源于天体计算来源于天体计算由若干观测值计算人一时刻星球的由若干观测值计算人一时刻星球的位置。现在，插值法在工程技术和数据处理有许多直接位置。现在，插值法在工程技术和数据处理有许多直接
2、应用，而且也是数值积分、数值微分的基础。应用，而且也是数值积分、数值微分的基础。2.1 插值概念与基础理论插值概念与基础理论2.1.1 插值问题的提法对于给定的函数表xx0 x1. xnY=f(x) y0y1. yn(其中其中 在在a,b上连续，上连续， x0，x1，xn 是是 a,b上的上的 n+1个互异的点个互异的点)，在某函数类，在某函数类 (x) 中求一个函数中求一个函数 (x) ，使，使( )yf x (xi)=yi , (i=0,1,2,n) (2)（1）并用函数并用函数 (x) 作为函数作为函数 y=f(x) 的近似函数，即的近似函数，即y= f(x) (x) , ( xa,b
3、) 这类问题称为这类问题称为插值问题插值问题。 a,b称为称为插值区间插值区间， x0 , x1, . , xn 称为称为插值节点插值节点，（，（2）称为）称为插值条件插值条件，插值条件，插值条件是选择近似函数的标准，满足此条件的近似函数是选择近似函数的标准，满足此条件的近似函数 (x) 称称为为插值函数插值函数， f(x) 称为称为被插值函数被插值函数。 函数类函数类 (x) 有多种取法，常用的有代数多项式、有多种取法，常用的有代数多项式、三角函数和有理函数。三角函数和有理函数。 最简单的插值函数是代数多项式最简单的插值函数是代数多项式，相应的插值问题称为多项式插值多项式插值。 最简单的插值
4、函数是代数多项式最简单的插值函数是代数多项式，相应的插值问题称为多项式插值多项式插值。 根据所给函数表（根据所给函数表（1），求一个次数不高于），求一个次数不高于n的多项的多项式式Pn(x)=a0+a1x+anxn, (3)使pn(xi)=yi, （ i= 0,1,2,，n） (4)满足插值条件（满足插值条件（4）的多项式（）的多项式（3），），称为函数称为函数y=f(x) 在在节点节点x0，x1，xn处的n次插值多项式。 2.2.2 多项式插值的理论基础多项式插值的理论基础 求)(xf的n次插值多项式( )nypx 的几何意义，就是)(xfy 上的若干个节点，作一条代数曲线( )nypx 来
5、近似代替曲线)(xfy 。如图所示。 通过曲线a0+a1x0+anx0n=y0 a0+a1x1+anx1n=y1 . (5)a0+a1xn+anxnn=yn 插值多项式的存在唯一性插值多项式的存在唯一性 由插值条件（由插值条件（4）知，插值多项式）知，插值多项式Pn(x)的系数的系数a0 ,a1, an满足下列线性方程组满足下列线性方程组 由于由于xi互异，所以互异，所以(6)右端不为零，从而方程组右端不为零，从而方程组(5)的的解解 a0 ,a1 ,an 存在且唯一。于是有存在且唯一。于是有 2n0002n1112nnnn1.1.D.1.xxxxxxxxx 0()ijj i nxx 而而ai
6、(i=0,1,2,n)的系数行列式是的系数行列式是Vandermonde行列式行列式(6) 定理定理1 满足插值条件（4）的n次次插值多项式是存在且唯一的。n 2.1.3误差估计误差估计 ( )nypx )(xfy 从前面的分析知道，用代数多项式从前面的分析知道，用代数多项式来近似代替曲线来近似代替曲线。除了在节点处没有误差外，除了在节点处没有误差外，在其它点上一般都有误差。若记在其它点上一般都有误差。若记 xpxfxRnn xpn xf的截断误差（插值余项）的截断误差（插值余项）则则 nRx就是用就是用代替代替定理定理2 2 设 , ,nfxC a b 0,x1,x2,xnx xpn xf则
7、对任意,bax为在n+1个节点，有余项 xnfxpxfxRnnnn1)1(!1)()()()(ba,其中其中 niinnxxxxxxxxx0101任意的 , ,xa b 1nfx 存在。上的n次插值多项式，（7）应当指出，余项表达式只有在应当指出，余项表达式只有在 f(x) 的高阶导数的高阶导数存在时才能应用。存在时才能应用。 在在 （a，b）内的具体位置通常不）内的具体位置通常不可能给出，如果我们可以求出可能给出，如果我们可以求出 ，那么插值多项式那么插值多项式Ln(x)逼近逼近f(x)的截断误差是的截断误差是(11) Mfnnbxax1)1)(max)()!1(11)(xnnMxRnn2.
8、2 插值多项式的求法插值多项式的求法 在前面讨论插值多项式的存在唯一性时，实际上已提供在前面讨论插值多项式的存在唯一性时，实际上已提供了它的一种求法，即通过求解线性方程组来确定其系数了它的一种求法，即通过求解线性方程组来确定其系数ai (i=0,1,2,n) 但是这种方法不仅计算量大，而且因不能获得简明的表但是这种方法不仅计算量大，而且因不能获得简明的表达式而给理论和应用研究带来不便。在这里我们学习两种达式而给理论和应用研究带来不便。在这里我们学习两种简便而实用的求答。简便而实用的求答。2.2.1 拉格朗日插值多项式拉格朗日插值多项式 在线性代数中知道，所有次数不超过在线性代数中知道，所有次数
9、不超过n次的多项式构次的多项式构成一个成一个n+1维线性空间。其基有各种不同的取法。因此维线性空间。其基有各种不同的取法。因此尽管满足条件（尽管满足条件（4）的）的n次插值多项式是唯一的，然而它次插值多项式是唯一的，然而它的表达式可以有多种不同的形式。如果取满足条件：的表达式可以有多种不同的形式。如果取满足条件：0,()1,kiikikl x 的一组n次多项式xlxlxlxln,210作为上述线性空间的基，则容易看出10010( )( )( )( )nnnk kklx ylx ylx yy lx n 是是一一个个次次数数不不超超过过 的的多多项项式式。且且满满足足插插值值条条件件（4 4）。因
10、此，由n+1个代数多项式 xlxlxlxln,210线性生成的多项式（10）就是满足插值条件的n次插值多项式。（10）（9）满足条件（9）的多项式称为n+1个节点的n次基本插值多项式（或n次基函数） xlxlxlxln,210 显然，求拉格朗日多项式的关键是求n次插值基函数。0,()1,kiikikl x 因此，可设 0111()().()().()kkkknlxA xxxxxxxxxx 因为 xlk为n次多项式，且n 011011()()()()( )()()()()nkkkkkkknkkxxxxxxxxlxxxxxxxxx 两种特殊的两种特殊的Lagrange插值多项式插值多项式1.线性插
11、值线性插值(两点插值两点插值) 最简单的插值是线性插值最简单的插值是线性插值(此时此时n=1), 这时插值问题就是这时插值问题就是求一次多项式求一次多项式P1(x)=a0+a1x 使它满足条件使它满足条件P1(x0)=y0 , P1(x1)=y1 ,这时这时1001( )xxlxxx 0110( )xxlxxx 于是线性插值多项式为于是线性插值多项式为011010110( )xxxxL xyyxxxx 即即100010( )()nyyLxyxxxx 它就是通过它就是通过M0(x0,y0)和和M1(x1,y1)两点的线段两点的线段.2.抛物插值抛物插值 线性插值仅仅用两个节点以上的信息，精确度较
12、差。为线性插值仅仅用两个节点以上的信息，精确度较差。为了提高精确度，我们进一步考察以下三点的插值问题了提高精确度，我们进一步考察以下三点的插值问题(n=2)：这时1200102()()( )()()xxxxlxxxxx 0211012()()( )()()xxxxlxxxxx 0122021()()( )()()xxxxlxxxxx 由此得到抛物插值多项式由此得到抛物插值多项式20 01 12 2( )( )( )( )L xy l xy l xy l x 抛物插值又称三点插值抛物插值又称三点插值.xyln 例例2 2 已知已知的函数表的函数表x 10 11 12 13 14y 2.3026
13、2.3979 2.4849 2.5649 2.6391 5 .11ln并估计误差。并估计误差。分别用拉格朗日线性和抛物线插值求分别用拉格朗日线性和抛物线插值求的近似值，的近似值，2.3 分段低次插值分段低次插值 插值的目的是数值逼近的一种手段，而数值逼近，为了得到一个数学问题的精确解或足够精确的解。那么，是否插值多项式的次数越高，越能够达到这个目的呢？现在，我们来讨论一下这个问题。 我们已经知道：f(x)在n+1个节点xi(i=0，1，2，n) 上的n次插值多项式Pn (x) 的余项 niinnxxnfxPxfxR011)()!()()()()()(设想当节点数增多时会出现什么情况。由插值余项
14、可知，当f(x)充分光滑时时，余项随n增大而减少，这说明可用增加节点的方法达到这个目的，那么实际是这样吗？1901年龙格(Runge) 给出一个例子: 定义在区间-1，1上，这是一个光滑函数，它的任意阶导数都存在，对它在-1，1上作等距节点插值时，插值多项式情况，见图:从图中，可见，在靠近-1或1时，余项会随n值增大而增大，如P12(0.96)=36!但f(0.96)=0.25 21()1fxx-5-4-3-2-1012345-1.5-1-0.500.511.52n=2n=4n=6n=8n=10f(x)=1/(1+x2) 从图中，还可发现，在0附近插值效果是好的，即余项较小，另一种现象是插值多
15、项式随节点增多而振动更多。 这种插值多项式当节点增加时反而不能更好地接近被插之数的现象，称为龙格现象龙格现象。 上述现象和定理，告诉我们用高次插值多项式是不妥当的，从数值计算上可解释为高次插值多项式的计算会带来舍入误差的增大，从而引起计算失真。那么如何提高插值精度呢？采用分段插值是一种办法。实践上作插值时一般只用一次、二次最多用三次插值多项式。分段线性插值的构造： 设f(x)是定义在a,b上的函数，在a,b上节点 a= x0 x1x2xn-1xn=b,的函数值为 y0 , y1 ,y2 ,yn-1 ,yn 。 (x)在每个子区间xi , xi+1(i=0,1,2,n-1)上是一次插值多项式;1
16、1111 iiiiiiiiiixxxxxxxyxxxxyx,)(这种分段低次插值称为分段线性插值分段线性插值.在几何上就是用折线段带代替曲线,故分段线性插值又称为折线插值折线插值.(,) ,0,1,kkxyin实际上是连接点的一条折线分段线性插值曲线图：-4-3-2-101234-1-0.8-0.6-0.4-0.200.20.40.60.81曲线的光滑性较差在节点处有尖点 但如果增加节点的数量减小步长,会改善插值效果分段二次插值分段二次插值即：选取跟节点x最近的三个节点xi-1,xi, xi+1进行二次插值，即在区间xi-1, xi+1，取：这种分段的低次插值叫分段二次插值，在几何上就是用分段
17、抛物线代替y=f(x)，故分段二次插值又和分段抛物插值。 11112iikikjijikjikxxxxyxLxf)()()()(什么是样条:是 指飞机或轮船等的制造过程中为描绘出光滑的外形曲线(放样)所用的工具样条本质上是一段一段的三次多项式拼合而成的曲线在拼接处,不仅函数是连续的,且一阶和二阶导数也是连续的1946年,Schoenberg将样条引入数学,即所谓的样条函数 2. 4 三次样条三次样条插值插值,)(,)(),(),()1(2baCxSbaxSxSxS 即上连续都在区间上都是三次多项式在每个小区间,)()2(1kkxxxS处的函数值为在节点如果函数nxxxxf,)()3(10njy
18、xfjj, 1 , 0,)( )S x而满足njyxSjj, 1 , 0,)(上的三次样条插值函数在为则称,)()(baxfxS-(1)定义1. 的一个分割为区间,10babxxxan:,)(上满足条件在区间如果函数baxS2.4.1、三次样条插值函数2( )s x由条件（ ），不妨将记为1 ,1,2,.,iiisxs x xx xin( )= ( ),32( )(1)iiiixa xb xc xdisiiiiabcd其中 ， ， ， 为待定系数，共4n个。1由条件（ ），有( )( )( )( )(1,2,.,1)(2)( )( )iiiiiixxxxinxxii+1ii+1ii+1ssss
19、ss3由条件（ ），有( )(0,1,2,., )(3)iiis xyin 2342n容易看出，（ ）和（ ）共有个方程，为确定s(x)的4n个待定参数，尚需再给出两个条件，即所谓边界条件。通常使用的边界条件有以下三类：000,nnnxfxfff第一类边界条件是s ()=s ()=为给定的值。00000( ) nnnnnf xx - xs xx - xs xs xs xs xs xs x第三类边界条件是周期条件。设 ( )是周期函数，不妨设以为一个周期，这时也应以为周期的周期函数，于是 s(x)在端点处满足条件：(+0)= (-0); (+0)= (-0); (+0)= (-0).000nnn
20、xfxfff第二类边界条件是s ()=s ()=， 为给定的值。0=0nxx当s ()=s ()时，样条函数在两端点不受力，呈自然状态，故称之为自然边界条件。4( )4( ).ns xns x 利用个条件求出三次样条函数的个待定常数，直接求解计算量很大，通常利用Matlab软件求例设f(x)为定义在0,3上的函数，有下列函数值表xi0123yi00.521.503()0.2,()1,0,3( ).fxfxs x 且试求区间上满足上述条件的三次养条插值函数解用Matlab求解s(x).程序为：x=0 1 2 3;y=0.2 0 0.5 2.0 1.5 -1;pp=csape(x,y,comple
21、te)breaks,coefs,npolys,ncoefs,dim=unmkpp(pp)( , ,)csape x y complete上述程序中函数是指定边界条件的样条插值函数，csape返回一个包含三次样条插值的pp形，或者说是分段多项式的结构。这个结构包含了计算用户希望的任何插值点数值的三次样条值需要的所有信息。字符串complete表示所给边界条件是第一类边界条件；若将complete换成second表示第二类边界条件；periodic表示第三类边界条件。在第一类边界条件和第二类边界条件时，边界条件值放在y的第一个分量和最后一个分量的位置上。这样y的分量的个数比x的分量的个数多。例如，
22、在例中，x=0 1 2 3y=0.2 0 0.5 2.0 1.5 -1周期边界条件时，无需指定边界条件值。在计算一个三次样条表达式的时候，必须将pp形中不同域提取出来进行计算，这个过程可以由函数unmkpp完成，该函数的使用方法为：breaks,coefs,npolys,ncoefs,dim=unmkpp(pp)其中输入变量pp是样条插值函数csape的输出变量，unmkpp的输出变量有个：breaks,coefs,npolys,ncoefs,dim。0123320111213,.( )()()():dim:iiiia a a as xa xxa xxaxxanpolysncoefsbreak
23、s：包含了插值节点；coefs：是一个矩阵，其第i行是第i个多项式的系数：即 中的系数。是多项式的个数；是每个多项式系数的个数；是样条的维数。运行结果如下：pp= form:ppbreaks:0 1 2 3coefs:3 4 doublepieces:3order:4dim:13232320.180.2 ,0,11.04(1)1.26(1)1.28(1)0.5, 1,22,3(2)1.86(2)0.68(2)2.0,xxxxs xxxxxxxxx因此所求的三次样条函数为0.48( )=0.68breaks= 0 1 2 3coefs= 0.2800 -0.1800 0.2000 0 -1.04
24、00 1.2600 1.2800 0.5000 0.6800 -1.8600 0.6800 2.0000npolys= 3ncoefs= 4dim= 1例已知函数值表xi1245yi1342( ).s x试求在区间,5 上满足上述函数表所给出的插值条件的三次样自然样条插值函数解用Matlab求解s(x).程序为：x=;y=;pp=csape(x,y,second)breaks,coefs,npolys,ncoefs,dim=unmkpp(pp)运行结果如下：pp= form:ppbreaks:1 2 4 5coefs:3 4 doublepieces:3order:4dim:1332321)2
25、.125(1) 1,1,20.125(2)0.375(2)1.75(2)3, 2,44,5(4)1.125(4)1.25(4)4,xxxs xxxxxxxxx因此所求的三次样条函数为-0.125( )=0.375breaks= 1 2 4 5coefs= -0.1250 0 2.1250 1.0000 -0.1250 -0.3750 1.7500 3.0000 0.3750 -1.1250 -1.2500 4.0000 在科学实验中，经常需要从一组实验数据( xi,yi)出发，求函数y=f(x)的一个近似表达式y=y(x)（通常称为经验公式）。从几何上看，就是通过给定m个数据点，求曲线y=f(
26、x)的一条近似曲线y=y(x)使这条曲线尽可能与所给的m个点相吻合。2.5 曲线拟合的最小二乘法曲线拟合的最小二乘法 例1 考察某种纤维的强度与其拉伸倍数的关系,下表是实际测定的24个纤维样品的强度与相应的拉伸倍数的记录:编 号 拉伸倍数 强 度编 号拉伸倍数 强 度11.91.41355.5221.3145.2532.11.81565.542.52.5166.36.452.72.8176.5662.72.5187.15.373.531986.583.52.72087944218.98.51043.52298114.54.2239.58.1124.63.524108.1iiyxiiyx 2.5
27、.1 最小二乘法的提法最小二乘法的提法 要求出拉伸强度和倍数的关系，插值法虽然在一定程度上，可以根据函数表求函数的近似表达式。但用来解决这里提出的问题存在明显缺陷： 1. 实验提供的数据带有误差，使用插值法会保留这些误差，从而失去原数据表示的规律。 2.实验数据往往很多，用插值法得到的近似表达式明显缺乏使用价值。 为了获得便于应用的经验公式，不用插值标准可能更合适，最小二乘法是解决这类问题的一种较好方法。1234567891012345678912345678910123456789系要关系应是线性关的主与拉伸倍数因此可以认为强度xy并且24个点大致分布在一条直线附近xxy10)(为待定参数其
28、中10,-(1)越接近越好样本点与所有的数据点我们希望),)()(10iiyxxxy必须找到一种度量标准来衡量什么曲线最接近所有数据点iiiyxy)(令一般使用mii0222在回归分析中称为残差miiiyxy02)(准偏离程度大小的度量标与数据点作为衡量),()(iiyxxy称为平方误差在回归分析中称为残差平方和从而确定(1)中的待定系数mii0222miiiyxy02)(残差的大小可以衡量近似函数的好坏。求出 使残差的平方和最小的方法称为曲线拟合的最小二乘法或最佳平方逼近。 , a b用最小二乘法求近似函数（经验公式）的基本步骤基本步骤如下： (1) 确定近似函数类,即确定近似函数 的形式。
29、 ( )yx 这并非单纯的数学问题,与其它各领域的专门知识有关.数学上,通常根据在坐标纸上所描点的情况来选择 ( )x 的形式 (2) 最小二乘解。 即求使残差的平方和最小的近似函数( )x 求系数 01,na aa 220111(,)nnmimiiiF a aayPx ,使 2.5.2 最小二乘法的求法最小二乘法的求法多项式拟合多项式拟合 ( ).f x( )mPx(,)(0,1, ),iixyin 设已知点求m次多项式 来拟合函数 mmmxaxaaxP10)(设达到最小。1020nmkjikiiikjFya x xa 1,1jm 11mnnkjjkiiik oiiaxy x 于是得线性方程
30、组即 当拟合函数是一元函数时，所对应的函数图形是平面曲线。这时，数据拟合问题的几何背景是寻求一条近似逼近给定离散点的曲线，故称为曲线拟合问题曲线拟合问题。 由于函数 11(,)mmF aaa达到最小,由高等数学知识有 111112111111112111111nnnmmmiiiiiinnnnmmimiiiiiiiinnnnmmmmmimiiiiiiiinaaxaxyaxaxaxx yaxaxaxx y -(3.1)*,(1,1)kkaakm*1*121( ).mmmmmPxa xa xa xa 方程组(3.1)称为法方程组法方程组.可以证明,该方程组有唯一解并且相应的函数 就是满足（3.1）的
31、最小二乘解最小二乘解。 21111231111123421111111221111.nnnmiiiiiinnnnmmiiiiiiiimnnnnmmiiiiiiiinnnnmmmmiiiiiiiinxxxaxxxxaaxxxxaxxxx11211.niiniiiniiinmiiiyx yx yx y 2.5.3 用用Matlab作最小二乘拟合作最小二乘拟合多项式在处的值可用下面程序计算：y=polyval(a,x)11 ( , ,);.,.,.mmmapolyfit x y mxxxa12m12m+1多项式拟合有现成程序：其中输入参数x,y是要拟合的数据，为长度自定义的数组，m为拟合多项式的次数，输出参数a为拟合多项式 y=aaa的系数 a=a aa例例2 设一发射源强度公式为 0,atII e 观测数据如下: 试用最小二乘法确定I与t的关系式。 解解 00,lnln.atII eIIat 将观测数据化为用Matlab编成如下：x=0.2:0.1:0.8;y=1.106 0.8671 0.5596 0.2827 0.0000 -0.3011 -0.5789;a=poly(x,y,1);(1)2.89, (2)1.73aa 运行程序所以 ln1.732.89 ,It即 005.64,2.89aIea 于是 2.895.64.tIe
