matlab插值法.ppt

第1章插值法,1.1插值法1.2Lagrange插值1.3Newton插值1.4Hermite插值1.5分段线性插值1.6三次样条插值1.7程序示例习题1,1.1插值法,插值问题的背景在生产和实验中，函数f(x)或者其表达式不便于计算,或者无表达式而只有函数在给定点的函数值(或其导数值)，此时我们希望建立一个简单的而便于计算的近似函数(x)，来逼近函数f(x)。常用的函数逼近方法有：插值法；最小二乘法（或称均方逼近）；一致逼近等。,插值法,插值法是函数逼近的重要方法之一，有着广泛的应用。简单地说，插值法就是用给定的(未知)函数f(x)的若干点上的函数值(或其导数值)来构造f(x)的近似函数(x)，要求(x)与f(x)在给定点的函数值相等。有很多种插值法，其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点，常用的插值还有Hermite插值，分段插值和样条插值。,函数可以未知，只需已知若干点上的值。,设f(x)为a，b上的函数，在互异点x0,x1,.,xn处的函数值分别为f(x0),f(x1),f(xn),构造一个简单函数(x)作为函数f(x)的近似表达式y=f(x)(x),使(xi)f(xi),i0,1,2,n(1.0)则称(x)为关于节点x0,x1,.,xn的插值函数；称x0,x1,.,xn为插值节点；称(xi,f(xi),i=1,2,n为插值点；f(x)称为被插值函数。(1.0)式称为插值条件。这类问题称为插值问题。构造出(x)，对f(x)在a,b上函数值的计算，就转化为(x)在对应点上的算。,插值法的定义,1.2Lagrange插值选用代数多项式作为插值函数。Lagrange插值就是选用节点上的函数值作为插值条件。,1.2.1线性插值给定两个点(x0，y0)，(x1，y1），x0x1，确定一个一次多项式插值函数，简称线性插值。待定系数法设L1(x)=a0+a1x,代入插值点当x0x1时，方程组的解存在唯一。,即插值条件：L1(xi)=f(xi)=yi,i=0,1,解之得，,因此，（1.1）式称为一次Lagrange插值。由求解过程知，用待定系数法，需要求解线性方程组，当已知节点较多时，即方程的未知数多，计算量较大，不便向高阶插值推广。,插值基函数法,分别构造两个节点上的一次函数，使其在本节点上的函数值为1，而在其他节点上的函数值为0。设l0(x),l1(x)分别为满足上述条件的一次函数，即或简单地记为对于过两个节点x0,x1的线性插值（1.1）式，令,显然，l0(x),l1(x)满足：线性插值函数可以写成节点上函数值的线性组合，即L1(x)=l0(x)y0+l1(x)y1称l0(x),l1(x)分别为x0,x1的插值基函数。,线性插值误差定理1设L1(x)为一次Lagrange插值函数,若f(x)一阶连续可导，f(x)在(a,b)上存在，则对任意给定的x(a,b),至少存在一点(a,b),使得证明因为L(xi)=f(xi),i=0,1,所以，R1(x0)=R1(x1)=0,即x0,x1为R1(x)的两个根。因此，可设R1(x)为,易知满足插值条件：L1(xi)=yi,i=0,1,可设R1(x)=k(x)(x-x0)(x-x1).,固定任一x，作辅助函数，令则(xi)=0,i=1,2,(x)=0,即(t)有3个零点x0,x1,x。假定，x0xx1,分别在x0，x和x，x1上应用洛尔（Rolle）定理，可知，(t)在每个区间上至少存在一个零点，1，2，使(1)=0，(2)=0(此即(t)有2个零点)。再利用洛尔定理知，(t)在1,2上至少有一个零点，使()=0。对(t)求2阶导数得，(t)=f(t)-2!k(x),因为()=0，所以，有k(x)=f()/2!。证毕。,1.2.2二次插值给定3个互异插值点(xi,f(xi),i=0,1,2，确定一个二次插值多项式函数，即抛物线插值（如图）。,待定系数法设L2(x)=a0+a1x+a2x2,代入3个插值条件：L2(xi)=f(xi),i=0,1,2，解线形方程组可得a0,a1,a2。,插值基函数法构造3个节点上2次插值基函数l0(x),l1(x),l2(x),使满足li(xj)=ij，i,j=0,1,2。因为l0(x)为2次插值基函数,且l0(x1)=l0(x2)=0,所以可设l0(x)=A(x-x1)(x-x2)。由条件：l0(x0)=1，得同理可得，二次Lagrange插值多项式为,容易验证满足插值条件,二次插值的误差,定理设L2(x)为二次Lagrange插值函数,若f(x)C3a，b，则任给x(a,b),至少存在一点=(x)(a,b),使提示：因为R2(x0)=R2(x1)=R2(x2)=0,可设作辅助函数易知，x0,x1,x2,x为(t)的4个零点，在4个点两两组成的区间上，应用Rolle定理，然后再反复应用Rolle定理即得证。,例1.1给定sin11=0.190809,sin12=0.207912,求线性插值，并计算sin1130和sin1030。,解x0=11,x1=12,y0=0.190809,y1=0.207912，sin1130L1(11.5)=0.199361,sin1030L1(10.5)=0.182258.由定理1知，误差为,准确值为：sin1130=0.199368sin1030=0.182236,例1.2给定sin11=0.190809,sin12=0.207912,sin13=0.224951,构造二次插值，并计算sin1130。,解x0=11,x1=12,x2=13,y0=0.190809,y1=0.207912,y2=0.224951，sin1130L2(11.5)=0.199369,sin1130=0.199368.,例1.3要制作三角函数sinx的值表，已知表值有四位小数，要求用线性插值引起的截断误差不超过表值的舍入误差，试确定其最大允许的步长。,解f(x)=sinx,设xi,xi为任意两个插值节点，最大允许步长记为h=hi=xixi，,1.2.3n次Lagrange插值多项式,已知n+1个互异插值节点(xi,f(xi),i=0,1,2,n,研究n次插值多项式的存在性及其表示形式。存在性设n次多项式为代入插值点，即插值条件：Pn(xi)=f(xi),i=0,1,2,n,得其范德蒙德(Vandermonde)行列式为：,所以，(1.6)的解存在唯一。解出(1.6)的解a0,a1,an，代入(1.6)1即得n次插值多项式Pn(x)。n次插值多项式的构造有上面的讨论知，用待定系数法要求解一个线性方程组，计算量大，实际中不可取。,插值基函数法,分别构造x0,x1,xn上的n次插值基函数l0(x),l1(x),ln(x),满足性质：即,先构造l0(x)。有上表知，x1,x2,xn为l0(x)的零点，设,由l0(x0)=1，得同理可设由li(xi)=1，得,于是，,所以我们得到n次Lagrange插值多项式：容易验证，Ln(xi)=f(xi),i=0,1,2,n.,例1.4已知插值点(-2.00,17.00),(0.00,1.00),(1.00,2.00),(2.00,17.00),求三次插值，并计算f(0.6)。,解先计算4个节点上的基函数：,三次Lagrange插值多项式为：f(0.6)L3(0.6)=-0.472.,n次插值多项式的误差定理2设Ln(x)是a,b上过插值点(xi,f(xi)的n次Lagrange插值多项式，节点xi两两互异，若f(x)Cn+1a,b,则xa,b,存在=(x)a,b,使得证明提示：因为Rn(xi)=0,i=0,1,2,n,令作辅助函数显然，(t)有n+2个零点x0，x1，xn，x，反复应用Rolle定理，由(n+1)()=0,定理得证。,Back,n次插值多项式的几点说明,若|f(n+1)(x)|M,xa,b,则由(1.9)得当f(x)为不高于n次的多项式时，因为Rn(x)=0,则Ln(x)=f(x).特别地，取f(x)=xk,k=0,1,2,n,则令k=0，可知n次Lagrange基函数li(x)满足,即f(x)为不超过n次的多项式时，插值多项式就是被插函数本身,实际计算中，如果f(x)的形式未知，也就无法求出f(n+1)(x)或估计出其较精确的上界。所以也就无法估计|Rn(x)|。我们采用事后估计，即给定n+2个节点，x0,x1,xn,xn+1,构造两个n次插值多项式Ln(x)和，用它们来估计|Rn(x)|。由定理2，得,Th.2,假定f(n+1)(x)在a,b内连续,且变化不大,即f(n+1)(1)f(n+1)(2).,(1.13)和(1.14)两式相除，得解之得，于是，得误差,事后估计：即用求出的插值多项式来估计误差。,n次插值多项式的算法描述,输入节点数n、插值点(xi,yi)(i=0,1,2,n)和要计算的函数点x。实现基函数li(x)和插值Ln(x):FORi=0,1,ntmp=1;FORj=0,1,i1,i+1,ntmp=tmp*(xxj)/(xixj);(实现li(x)fx=fx+tmp*yi;(实现Ln(x)输出fx。,2006年9月22日,1.3Newton插值,Lagrange插值的优缺点：优点：形式整齐、规范，理论上保证插值的存在唯一性。缺点：计算量大、不具有承袭性。1.3.1差商及其性质一阶差商：f(x)关于点x0，x1的一阶差商记为fx0,x1，二阶差商：f(x)关于点x0，x1，x2的二阶差商记为fx0,x1,x2,一般地，k阶差商fx0,x1,xn定义为：,差商的性质性质1k阶差商fx0,x1,xk可表成节点上函数值f(x0),f(x1),f(xk)的线性组合，即例如，k=2时，(1.17),性质2各阶差商具有对称性,即改变差商中节点的次序不会改变差商的值。设i0,i1,ik为0,1,k的任一排列,则,由性质1知，任意改变节点的次序，只改变(1.17)式右端求和的次序，故其值不变。例如,由定义知，性质3若f(x)为n次多项式，则一阶差商fx,xi为n1次多项式。由定义令x=xi,则分子为0,说明分子中含有因子xxi,与分母约去。,性质4若f(x)在a,b存在n+1阶导数，xia,b,i=0,1,n,固定xa,b,则n+1阶差商与导数存在如下关系：,差商的计算由差商的定义一阶差商是由节点上函数值定义的，二阶差商是由一阶差商定义的，依此构造差商表：,Return,例1.5计算(2,17),(0,1),(1,2),(2,19)的一至三阶差商。,解由表易知，fx0,x1=f2,0=8fx0,x1,x2=f2,0,1=(81)/(21)=3fx0,x1,x2,x3=f2,0,1,2=(38)/(22)=5/4,例1.7,由差商与导数的关系（性质4）,例1.6对f(x)=x7x4+3x+1，求f20,21，fx,20,21,26和fx,20,21,27。解显然，f(7)(x)=7!,f(8)(x)=0,由性质4得,2006.9.22,1.3.2Newton插值,线性插值给定两个插值点(x0,f(x0),(x1,f(x1),x0x1,设N1(x)=a0+a1(xx0)直线的点斜式代入插值点得，于是得线性Newton插值公式,由插值的唯一性知，L1(x)与N1(x)为同一多项式，只是表达形式不同而已。,二次Newton插值,给定三个互异插值点(xi,f(xi),设代入插值条件：N2(xi)=f(xi),i=0,1,2,得二次Newton插值公式为,n次Newton插值公式,给定n+1个插值点(xi,f(xi),i=0,1,2,n,xi互异，类似地，有二阶至n阶差商的定义得(xx0)(xx0)(xx1)(xx0)(xxn1)上述所有n+1个等式相加，得,其中，插值误差为：,n次Newton插值公式,容易验证，Newton插值满足插值条件：Nn(xi)=f(xi),i=0,1,2,n.,关于Lagrange插值和Newton插值的几点说明由插值的唯一性质，Ln(x)=Nn(x)。因此，他们的误差也相同，即当f(x)Cn+1a,b时，有故得差商的性质4,2.牛顿插值的误差不要求函数的高阶导数存在，所以更具有一般性。它对f(x)是由离散点给出的函数情形或f(x)的导数不存在的情形均适用。,3.引入记号:fx0=f(x0),t0(x)=1,t1(x)=xx0,t2(x)=(xx0)(xx1),tn(x)=(xx0)(xx1)(xxn1)，于是n次Newton插值公式可表为称t0(x),t1(x),t2(x),tn(x)为Newton插值的基函数，而且满,足如下关系：ti(x)=ti1(x)(xxi1),i=1,2,n；ti(xj)=0，ji，ti(xj)0，ji。,4.Newton插值具有承袭性质，即5.Newton插值公式的计算量乘：1+2+(n1)+n=n(n+1)/2除：n+(n1)+2+1=n(n+1)/2,第i个节点以后非零,例1.7给定四个插值点(2,17),(0,1),(1,2),(2,19),计算N2(0.9),N3(0.9)。,解x0=2，x1=0，x2=1，x3=2，由例1.5知，fx0,x1=8，fx0,x1,x2=3，fx0,x1,x2,x3=5/4，所以，N2(0.9)=178(0.9+2)+3(0.9+2)0.9=1.63;N3(0.9)=N2(0.9)+1.250.9(0.9+2)(0.91)=1.30375.,例1.5,Newton插值的算法,用牛顿插值公式，首先要计算各阶差商值，然后再计算插值。牛顿插值由承袭性质容易实现，关键是实现差商。1.输入初始数据:节点数n,插值点序列x，f(xi),i=0,1,2,n,及要计算的函数点u；2.形成差商表gk,k=1,2,n;（gk=fx0,xk）3.形成插值基函数及插值t=1,newton=f(x0)对i=1,2,nt=(uxi1)t;(由tk=(uxk1)tk1形成(ux0)(uxi1)newton=newton+tgi4.输出牛顿插值Nn(u)=newton。,牛顿插值公式中只用到差商表中对角线上的差商，即fx0,fx0,x1,fx0,x1,x2,fx0,x1,xn.可以分别用一维数组gi来存放这些差商，i=0,1,2,n.,形成差商具体步骤：(1)对gi初始化，即gi=f(i),i=0,1,2,n.(2)除了g0(即f(x0)以外，其余函数值在计算一阶差商后不再使用，因此可用来存放一阶差商，即gj=fxj1,xj,j=n,n1,2,1(3)类似地，计算二阶差商后，除g1=fx0,x1外,其余一阶差商也不再使用,可用gj存放二阶差商fxj2,xj1,xj,即gj=fxj2,xj1,xj,j=n,n1,2，(4)最后，gn=fx0,x1,xn.,差商表,程序1,计算差商算法：fori=0tongi=f(xi)fork=1tonforj=ntokgj=(gjgj1)/(xjxjk),1.4埃尔米特(Hermite)插值,Hermite插值描述：设f(x)具有一阶连续导数，已知节点上的函数值和导数值，即(xi,f(xi),(xi,f(xi),i=0,1,2,n,若存在2n+1次多项式H2n+1(x)满足则称H2n+1(x)为f(x)关于节点xi(i=0,1,2,n)的Hermite插值多项式。记f(xi)=yi,f(xi)=mi,i=0,1,2,n.,三次Hermite插值的构造,存在性给定f(xi)=yi,f(xi)=mi,i=0,1.设代入插值条件:H3(xi)=f(xi),H3(xi)=f(xi),i=0,1,得,其解存在唯一,解出a0,a1,a2,a3,代入即得H3(x).,基函数法,每个节点上对应两套基函数：x0:h0(x),g0(x);x1:h1(x),g1(x),满足或简记为,先构造h0(x),设h0(x)=(a+bx)(xx1)2,为方便计算，可设,由h0(x0)=1,得a=1；由所以，同理设,h0(x1)=h0(x1)=0,g0(x0)=g0(x1)=0,g0(x1)=0,由g0(x0)=1,得a=1。所以,注：我们知道，过x0,x1两点的Lagrange插值基函数为显然，于是，三次Hermite插值的基函数可表为,到2n+1次,三次Hermite插值多项式为容易验证，当f(x)C4a,b时,三次Hermite插值的误差为提示：设作辅助函数固定xa,b，则(t)有三个零点x0,x1,x,且x0,x1为二重零点。反复应用Rolle定理可证。,2005年9月30日,高次Hermite插值的构造插值基函数法,给定n+1个节点x0,x1,xn上的函数值f(xi)和导数值f(xi)，可以构造2n+1次Hermite插值多项式满足其中，hi(x),gi(x)分别为对应于函数值和导数值的不高于2n+1次插值基函数，它们满足,完全仿照三次Hermite插值基函数的求法，可得,容易验证，Hermite插值的误差(定理)：当f(x)C2n+2a,b时,则存在(a,b),使,提示：对li(x)取对数ln，并注意li(xi)=1,见三次基,例1.8给定f(1)=0,f(1)=4,f(1)=2,f(1)=0,求H3(x),并计算f(0.5).,解x0=1,x1=1,f(0.5)H3(0.5)=3.5625.,例1.9给定f(0)=1,f(1)=2,f(0)=2,构造二次插值函数。,解(1)公式法设f(1)=m1，有三次Hermite插值公式得，令m1=0，得到二次Hermite插值函数P2(x)=x2+2x+1.,(2)插值基函数法,设节点x0上的二次基函数为t0(x),g0(x),节点x1上的二次基函数为t1(x),它们满足设由t0(x0)=1,即b(01)=1,得b=1。由t0(x0)=0,即a(01)+a0+b=0,得a=1。所以同理可设，分别由条件,(3)扩展牛顿法,写成差商表的形式，将带导数的节点X0及其上的函数值重复一遍，无导数的节点X1不重复，即xf(x)f(x)x001x1012X1x21211,扩展牛顿法用牛顿差商表构造Hermite插值,给定插值点(xi,f(xi),f(xi),i=1,2,n,重新定义插值节点序列：z2i=z2i+1=xi,i=0,1,2,n,即函数值取相应节点上的函数值,即f(z2i)=f(z2i+1)=f(xi),i=0,1,2,n；对于一阶差商，取,以偶数节点开始,以奇数节点开始,二阶以后的各阶差商，直接按差商公式计算。由此得到差商型Hermite插值公式：,差商表：zf(z)fzi,zj,例1.10已知计算f(1.36)。,解xf(x)f(x)1.20.61.20.60.51.40.91.551.40.90.74451.61.11.01.513145,1.5分段线性插值,n次Lagrange插值多项式的误差：插值多项式与被插函数的逼近程度同分点的数目和位置有关。一般地，分点越多，逼近程度越好，但也有例外。例如将1,110等分，步长h=2/10=0.2,取节点xi=1+0.2i,i=0,1,2,10。以(xi,f(xi)为插值点，构造L10(x)：,图示,Runge现象插值多项式在插值区间上发生剧烈的震荡。它揭示了高次插值多项式存在的缺陷。产生的原因误差有截断误差和舍入误差两部分组成，而在插值的计算过程中，舍入误差可能会扩散或放大。,返回,分段线性插值,给定n+1个插值点：(xi,f(xi),i=0,1,2,n,在每个小区间xi,xi+1上作线性插值，节点xi,xi+1上的基函数分别为：显然满足分段线性插值为：xi,xi+1上,区间a,b上的线性插值p(x)就是将每个小区间xi,xi+1上的线性插值pi(x)连接起来,p(x)为xi,xi+1上不高于一次的多项式。即,p(x)的图形是一条以(xi,f(xi)为折点的折线。,用分段线性插值逼近上述例子的效果，取n=10。,Runge现象,例1.11已知计算f(1.2),f(3.3).,解,1.6三次样条插值函数,分段线性插值：已知节点上的函数值，插值函数整体续。三次Hermite插值：已知节上的函数值和导数值，插值函数具有一阶连续的导数。为了得到光滑度更高的插值函数，引入样条插值函数。“样条”名词来源于工程中船体和汽车等的外形设计：给出外形曲线上的一组离散点(样点)，如(xi,yi)，i=0,1,2,n,将有弹性的细长木条或钢条(样条)在样点上固定，使其在其它地方自由弯曲，这样样条所表示的曲线，称为样条曲线(函数)。,在数学上，它表现为近似于一条分段的三次多项式，它要求在节点处具有一阶和二阶连续导数。,定义给定a,b上n+1个节点：a=x0x1xn=b及节点上的函数值f(xi)=yi,i=0,1,2,n,若S(x)满足：(1)S(xi)=yi,i=0,1,n;(2)S(x)在xi,xi+1上至多是一个三次多项式；(3)S(x)C2a,b.则称S(x)为f(x)关于节点x0,x1,xn的三次样条插值函数，称xi为样条节点。,在xi,xi+1上构造一个三次多项式,共有4n个未知量，因此需要4n个条件。由定义中的(1)知，有n+1个条件；由定义中的(3)知，有3n3个条件，即再附加2个边界条件，就可得到4n个条件。三次分段样条插值函数可唯一确定。M关系式：用节点处的二阶导数表示样条插值函数。m关系式：用节点处的一阶导数表示样条插值函数。,1.6.1M关系式,引入记号：S(x)为xi,xi+1上的三次多项式,S(x)为一次函数,我们用节点上的二阶导数值Mi表示线性函数。设其中，对(1.29)积分两次,得,将S(xi)=yi,S(xi+1)=yi+1,代入(1.30)得到二元一次方程组,解得,将C，D代入(1.30)式，得到xi,xi+1上的三次样条插值函数,在xi点，左导数=右导数,Back例,令(1.32)式为n+1个未知数M0,M1,Mn的n-1个方程组，补充两个边界条件后，方程组(1.32)就有唯一解。（2005-10-9）,返回例1.12,程序2,下面给出三种边界条件：(1)给定M0,Mn的值，可以得到n-1个未知量的n-1个方程的方程组,对角占优的三对角带状矩阵,M0=0,Mn=0时，称为自然边界条件,back,(2)给定条件：S(x0)=m0,S(xn)=mn,分别由x0,x1,xn-1,xn上的三次样条插值函数此二式和(1.32)组成n+1个方程的方程组：,（3）若f(x)为以xn-x0为周期的周期函数，即y0=yn，则要求S(x)也为周期函数，即S(x0)=S(xn),于是有即m0=mn,M0=Mn,此时(1.32)变为n个未知量、n个方程的方程组。,即为关于M0,M1,Mn的等式,例1.12已知离散点:(1.1,0.4000),(1.2,0.8000),(1.4,1.6500),(1.5,1.8000),取自然边界条件M0=Mn=0,构造三次样条插值函数，并计算f(1.25).,解n=3.h0=x1-x0=0.1,h1=0.2,h2=0.1,