插值法及其matlab实现教案.pptx

会计学1插值法及其插值法及其matlab实现实现 描述事物之间的数量关系：函数。有两种情况：一是表格表格形式一组离散的数据离散的数据来表示函数关系；另一种是函数虽然有明显的有明显的表达式表达式，但很复杂，但很复杂，不便于研究和使用。从实际需要出发：对于计算结果允许有一定的误差，可以把函数关系用一个简单的便于计算和处理的近似表达式近似表达式来代替，从而使问题得到简化。一般地，构造某种简单函数代替原来函数。0 引言引言第二章第二章 插值插值（Interpolation）法法插值法就是一种基本方法插值法就是一种基本方法第1页/共160页当精确函数当精确函数 y=f(x)非常复杂或未知时，在一系列节非常复杂或未知时，在一系列节点点 x0 xn 处测得函数值处测得函数值 y0=f(x0),yn=f(xn)，由，由此构造一个简单易算的近似函数此构造一个简单易算的近似函数 g(x)f(x)，满足条，满足条件件g(xi)=f(xi)(i=0,n)。这里的。这里的 g(x)称为称为f(x)的的插值函数插值函数。x0 x1x2x3x4xg(x)f(x)第2页/共160页 根据实际需要，可以用各种不同的函数来近似原来的函数。最常用的插值函数是最常用的插值函数是?多项式：多项式：代数多项式最简单，计算其值只需用到加、减乘运算，且积分和微分都很方便；所以常用它来近似表示表格函数(或复杂函数),这样的插值方法叫做代数插值法代数插值法，简称插值法。第3页/共160页1 拉格朗日多项式拉格朗日多项式 niyxPiin,.,0,)(=求求 n 次多项式次多项式 使得使得条件：条件：无重合节点，即无重合节点，即n=1已知已知 x0,x1;y0,y1，求，求使得使得111001)(,)(yxPyxP=可见可见 P1(x)是过是过(x0,y0)和和(x1,y1)两点的直线。两点的直线。)(1xP101xxxx 010 xxxx =y0+y11.1线性插值线性插值两点两点式式)()(0010101xxxxyyyxP +=点斜点斜式式)(001010 xxxxxxy +=()ff第4页/共160页1.2 二次插值二次插值n=2已知已知 x0,x1,x2;y0,y1,y2,求求使得使得002,)(yxP112)(yxP=222)(yxP=,为求为求P2(x),将三点代入其表达式将三点代入其表达式,即可得到三个方程式即可得到三个方程式,从而联立方程组解出系数从而联立方程组解出系数a0,a1,a2即可即可:方程组的方程组的解是否存解是否存在在?若存在解若存在解,是否是否唯一唯一?!当当 x0,x1,x2互异时互异时,方程组的解存在且唯一方程组的解存在且唯一.注：注：显然有显然有,求求n 次插值时次插值时,由由n+1个点可有个点可有n+1个方程个方程,联立方程组即可求出插值多项式的联立方程组即可求出插值多项式的n+1个系数个系数.然而然而,方程组的求解也并不是一件容易的事方程组的求解也并不是一件容易的事。待定系数法待定系数法第5页/共160页 对于线性插值的对于线性插值的两种形式两种形式解解进行适当的分析进行适当的分析,从中寻从中寻求规律而得到启发求规律而得到启发,就有了所谓的就有了所谓的拉格朗日插值法拉格朗日插值法(公式公式)和和牛顿插值牛顿插值(公式公式).).我们先来看看如何得到二次二次拉格朗日插值公式。拉格朗日插值公式。第6页/共160页 首先,线性插值的两点式可看作是两个特殊的一次式可看作是两个特殊的一次式的一种线性组合的一种线性组合.101xxxx 010 xxxx )(1xP=y0+y1=10)(iiiyxl两点两点式式l0(x)l1(x)实质上实质上（）和）和（）即是满足函数表）即是满足函数表 的一次插值多项式的一次插值多项式,称称l0(x)和和l1(x)为以为以x0,x1为节点的基本插值多项式，也称为为节点的基本插值多项式，也称为线性插值的线性插值的插值基函数插值基函数。于是，线性插值即是用于是，线性插值即是用基函数的线性组合基函数的线性组合来构造的来构造的.基函数法基函数法称为称为拉氏基函拉氏基函数数 ，满足，满足 li(xj)=ij 显然有显然有l0(x)+l1(x)1.这里，这里，l0(x)和和l1(x)具有如下性质：具有如下性质：l0(x0)=1,l0(x1)=0,l1(x0)=0,l1(x1)=1,第7页/共160页 由此启发，我们希望二次插值也能由一些二次插值基函数来线性组合由此启发，我们希望二次插值也能由一些二次插值基函数来线性组合:这时，这时，l0(x),l1(x),l2(x)都是二次多项式，且应满足都是二次多项式，且应满足满足满足(2.1)式的式的 l i(x)是否存在是否存在?若存在，具有什么形式呢若存在，具有什么形式呢?（2.1）第8页/共160页同理可得同理可得 l1(x)1(x x0)(x x2),l2(x)2(x x0)(x x1),1(x1x0)(x1x2)12(x2x0)(x2x1)1此即此即二次二次拉格朗日插值公式拉格朗日插值公式,其中其中,l0(x),l1(x),l2(x)是满足是满足(2.1)的特殊的特殊(基本基本)二次插值多项式二次插值多项式;称为称为二次插值基函数二次插值基函数.P2(x)=y0+y1+y2(x-x0)(x-x2)(x1-x0)(x1-x2)(x-x1)(x-x2)(x0-x1)(x0-x2)(x-x0)(x-x1)(x2-x0)(x2-x1)先考虑先考虑 l0(x)。因。因 l0(x)是以是以 x1,x2 为零点的二次多项式为零点的二次多项式,所以它可写成所以它可写成 l0(x)0(x x1)(x x2),其中其中 0 是待定系数。是待定系数。又因为又因为 l0(x0)=1，所以，所以 0(x0 x1)(x0 x2)1，则可有，则可有0(x0 x1)(x0 x2)1 l0(x)0(x x1)(x x2),第9页/共160页n 1希望找到希望找到li(x)，i=0,n 使得使得 li(xj)=ij；然后令；然后令=niiinyxlxP0)()(，则显然有，则显然有Pn(xi)=yi。li(x)每个每个 li 有有 n 个根个根 x0 xi xn=njj i jiniiixxCxxxxxxCxl00)().().()(=j i jiiiixxCxl)(11)(拉格朗日拉格朗日 多项式多项式与与 有关，而与有关，而与 无关无关节点节点f1.3 n 次插值次插值第10页/共160页定理定理 (唯一性唯一性)满足满足 的的 n 阶插值阶插值多项式是唯一存在的。多项式是唯一存在的。证明：证明：(存在性存在性可利用可利用Vandermonde 行列式行列式论证论证)反证：若不唯一，则除了反证：若不唯一，则除了Ln(x)外还有另一外还有另一 n 阶多项阶多项式式 Pn(x)满足满足 Pn(xi)=yi。考察考察 则则 Qn 的阶数的阶数 n而而 Qn 有有 个不同的根个不同的根n+1x0 xn注：注：若不将多项式次数限制为若不将多项式次数限制为 n，则插值多项式，则插值多项式不唯一不唯一。例如例如 也是一个插值多项式，也是一个插值多项式，其中其中 可以是任意多项式。可以是任意多项式。第11页/共160页设节点设节点在在a,b内存在内存在,考察截断误差考察截断误差，且，且 f 满足条件满足条件 ,Rolles Theorem:若若 充分光滑，充分光滑，则，则存在存在 使得使得 。推广：推广：若若使使得得使使得得存存在在使使得得Rn(x)至少有至少有 个根个根n+1=niinxxxKxR0)()()(任意固定任意固定 x xi (i=0,n),考察考察=niixtxKtRnt0)()()()(t)有有 n+2 个不同的根个不同的根 x0 xn x!)1()()()1(+-+nxKRxnn 注意这里是对注意这里是对 t 求导求导=+!)1)()()()1()1(nxKLfxnnxn !)1()()()1(+=+nfxKxn 1.4 插值余项插值余项(Remainder)第12页/共160页注：注：通常不能确定通常不能确定 x,而是估计而是估计 ,x(a,b)将将 作为误差估计上限。作为误差估计上限。当当 f(x)为任一个次数为任一个次数 n 的的多项式多项式时，时，,可知可知 ，即插值多项式对于次数，即插值多项式对于次数 n 的的多项式多项式是是精确精确的。的。第13页/共160页例例1 求经过求经过A(0,1)，B(1,2)，C(2,3)三个插值点的插值多项式三个插值点的插值多项式.解：解：三个插值节点及对应的函数值为三个插值节点及对应的函数值为由抛物插值公式得由抛物插值公式得第14页/共160页例例2：已已知知分别利用分别利用 sin x 的的1次、次、2次次 Lagrange 插值计算插值计算 sin 50 并估计误差。并估计误差。解：解：n=1分别利用分别利用x0,x1 以及以及 x1,x2 计算计算利用利用这里这里而而sin 50 =0.7660444)185(50sin10 p pL0.77614外推外推 (extrapolation)的实际误差的实际误差 0.010010.01001利用利用sin 50 0.76008,内插内插 (interpolation)的实际误差的实际误差 0.005960.00596内插通常优于外推。选择要计算的内插通常优于外推。选择要计算的 x 所在的区间的端点，插值效果较好。所在的区间的端点，插值效果较好。第15页/共160页n=2)185(50sin20 p pL0.76543sin 50 =0.76604442次插值的实际误差次插值的实际误差 0.000610.00061高次插值通常优于低次插值高次插值通常优于低次插值,但绝对不是次数越高就越好但绝对不是次数越高就越好第16页/共160页 例例3 考虑下述的插值法问题：求二次多项式考虑下述的插值法问题：求二次多项式P(x)，满足，满足P(x0)=y0,其中其中 是已给的数据并给出使这一问题的解存在且唯一的条件是已给的数据并给出使这一问题的解存在且唯一的条件.解解：设：设 则则 由已知条件有由已知条件有即即 所以所以故原问题的唯一可解性就归结为上述方程组的唯一可解性而后者唯一可解的充要条件为故原问题的唯一可解性就归结为上述方程组的唯一可解性而后者唯一可解的充要条件为这就是这就是P（x）存在且唯一的条件。）存在且唯一的条件。第17页/共160页 回顾回顾 拉格朗日插值公式拉格朗日插值公式 1.5 拉格朗日插值公式的优缺点拉格朗日插值公式的优缺点第18页/共160页n 1希望找到希望找到li(x)，i=0,n 使得使得 li(xj)=ij；然后令；然后令=niiinyxlxP0)()(，则显然有，则显然有Pn(xi)=yi。li(x)每个每个 li 有有 n 个根个根 x0 xi xn=njj i jiniiixxCxxxxxxCxl00)().().()(=j i jiiiixxCxl)(11)(拉格朗日拉格朗日 多项式多项式与与 有关，而与有关，而与 无关无关节点节点f拉格朗日插值公式拉格朗日插值公式 第19页/共160页 Lagrange插值插值公式公式(利用利用插值基函数插值基函数很容易得到很容易得到):含义直观含义直观,结构紧凑结构紧凑,在理论分析中非常方便在理论分析中非常方便;计算机上计算机上实现实现也很也很容易容易.也有一些也有一些缺点缺点：一是一是计算量大计算量大，这是显然的；另外，还有一个更严重的缺点，当插值节点增加时，这是显然的；另外，还有一个更严重的缺点，当插值节点增加时，全部插值基函数全部插值基函数均要随之均要随之变化变化，整个计算工作必须从头开始：不仅原来的每一项都要改变，还要，整个计算工作必须从头开始：不仅原来的每一项都要改变，还要增加一项增加一项计算。计算。为克服上述两个缺点为克服上述两个缺点,努力：把插值多项式变形为努力：把插值多项式变形为便于计算便于计算的形式。的形式。希望：希望：计算改变的过程中计算改变的过程中,尽可能能利用已有的计算结果尽可能能利用已有的计算结果.下面我们将看到下面我们将看到,这是可能的。我们可以有具有这是可能的。我们可以有具有“承袭性承袭性”的所谓牛顿公式。的所谓牛顿公式。第20页/共160页*1 拉格朗日插值的拉格朗日插值的Matlab实现实现常用的几个函数常用的几个函数(一一)POLY 函数函数调用格式调用格式：Y=poly(V)(二二)POLYVAL 函数函数调用格式调用格式：Y=polyval(p,x)返回矩阵特征多项式的系数。返回n次多项式p在x处的值。输入变量p=p0 p1 p2pn是一个长度为n+1的向量，其元素为按降排列的多项式系数。如输入程序如输入程序Y=poly(3)；则返回；则返回1 -3，即其特征多项式为，即其特征多项式为x-3如如对多项式对多项式 ，计算在计算在x=5,7,9 的值。的值。第21页/共160页(三三)POLY2SYM 函数函数调用格式一调用格式一：poly2sym(C)调用格式二调用格式二：f1=poly2sym(C,V)或 f2=poly2sym(C,sym(V)输入程序 p=3 2 1;x=5,7,9;polyval(p,x)结果为ans=86 162 262把多项式的系数向量转化为符号多项式如输入程序如输入程序 c=1-3;y=poly2sym(c)；则返回；则返回y=x-3如输入程序如输入程序 c=1-3;y=poly2sym(c,s)；或者或者 c=1-3;y=poly2sym(c,sym(s)；则返回则返回y=s-3第22页/共160页(四四)CONV 函数函数调用格式调用格式：C=conv(A,B)返回两个多项式乘积的多项式系数，即卷积和如输入程序如输入程序 A=1-3;B=1-2;C=conv(A,B)则返回则返回C=1 -5 6A、B为两个多项式的系数(五五)DECONV 函数函数调用格式调用格式：Q,R=deconv(B,A)conv函数的逆函数，返回两个多项式相除的多项式系数及其余数，即反卷积和若输入程序若输入程序 C=1-5 7;A=1-3;B,R=deconv(C,A)则输出则输出B=1 -2 R=0 0 1如输入程序如输入程序 C=1-5 6;A=1-3;B,R=deconv(C,A)则输出则输出B=1 -2 R=0 0 0第23页/共160页(六六)roots(poly(1:n)命令命令调用格式：调用格式：roots(poly(1:n)(七七)det(a*eye(size(A)-A)命令命令调用格式：调用格式：b=det(a*eye(size(A)-A)如输入程序如输入程序 C=1-3 2;roots(C)输出结果为输出结果为ans=2.0000 1.0000数值解数值解输入输入 poly(1:2)输出输出ans=1 -3 2输入输入 roots(poly(1:2)输出输出ans=2 11到到n的等差数列，其中的等差数列，其中n roots(poly(1:4)输出输出ans=4.0000 3.0000 2.0000 1.0000若常数a为矩阵A的特征值，则返回0返回n到1的数值第24页/共160页例例*1-1已知函数f(x)在1,3上具有二阶连续导数，且满足条件f(1)=1，f(3)=2，求线性插值多项式和函数值 f(1.5)，并估计其误差。解解：输入程序 X=1,3;Y=1,2;l01=poly(X(2)/(X(1)-X(2),l11=poly(X(1)/(X(2)-X(1),l0=poly2sym(l01),l1=poly2sym(l11),P=l01*Y(1)+l11*Y(2),L=poly2sym(P),x=1.5;Y=polyval(P,x)运行后输出基函数l0和l1及其插值多项式的系数向量P、插值多项式L和插值Y为l0=l1=L=Y=-1/2*x+3/2 1/2*x-1/2 1/2*x+1/2 1.2500一、线性插值的一、线性插值的Matlab程序程序逗号会将本句程序的结果输出到命令窗口，分号则不会输出第25页/共160页输入程序 M=5;R1=M*abs(x-X(1)*(x-X(2)/2运行后输出误差限为 R1=1.8750例例*1-2求函数f(x)=e-x在0,1上线性插值多项式，并估计其误差.解：解：输入程序 X=0,1;Y=exp(-X),l01=poly(X(2)/(X(1)-X(2),l11=poly(X(1)/(X(2)-X(1),l0=poly2sym(l01),l1=poly2sym(l11),P=l01*Y(1)+l11*Y(2),L=poly2sym(P),运行后输出基函数l0和l1及其插值多项式的系数向量P和插值多项式L为l0=-x+1 l1=x P=-0.6321 1.0000第26页/共160页输入程序 M=1;x=0:0.001:1;R1=M*max(abs(x-X(1).*(x-X(2)./2运行后输出误差限为 R1=0.1250.点乘和点除表示向量对应相乘和对应相除二、抛物线插值的二、抛物线插值的Matlab程序程序例例*1-3求将区间0,/2分成n等份（n=1,2），用y=f(x)=cos(x)产生n+1个节点，分别作线性插值函数和抛物线插值函数，并用它们分别计算cos(/6)(取四位有效数字)，并估计其误差.解：解：n=1输入程序 X=0,pi/2;Y=cos(X),l01=poly(X(2)/(X(1)-X(2),l11=poly(X(1)/(X(2)-X(1),l0=poly2sym(l01),l1=poly2sym(l11),P=l01*Y(1)+l11*Y(2),L=poly2sym(P),x=pi/6;Y=polyval(P,x)第27页/共160页运行后输出基函数l0和l1及其插值多项式的系数向量P、插值多项式和插值为P=-0.6366 1.0000Y=0.6667输入程序 M=1;x=pi/6;R1=M*abs(x-X(1)*(x-X(2)/2运行后输出误差限为R1=0.2742.n=2输入程序 X=0:pi/4:pi/2;Y=cos(X),l01=conv(poly(X(2),poly(X(3)/(X(1)-X(2)*(X(1)-X(3),l11=conv(poly(X(1),poly(X(3)/(X(2)-X(1)*(X(2)-X(3),l21=conv(poly(X(1),poly(X(2)/(X(3)-X(1)*(X(3)-X(2),第28页/共160页l0=poly2sym(l01),l1=poly2sym(l11),l2=poly2sym(l21),P=l01*Y(1)+l11*Y(2)+l21*Y(3),L=poly2sym(P),x=pi/6;Y=polyval(P,x)运行后输出基函数l01、l11和l21及其插值多项式的系数向量P、插值多项式L和插值Y为P=-0.3357 -0.1092 1.0000Y=0.8508第29页/共160页输入程序 M=1;x=pi/6;R2=M*abs(x-X(1)*(x-X(2)*(x-X(3)/6运行后输出误差限为R2=0.0239.cos(/6)=0.8660254三、三、n次拉格朗日插值的次拉格朗日插值的Matlab程序程序例例*1-4给出节点数据f(-2.00)=17.00，f(0.00)=1.00，f(1.00)=2.00，f(2.00)=17.00，作三次拉格朗日插值多项式计算f(0.6)，并估计其误差。解解：输入程序 X=-2,0,1,2;Y=17,1,2,17;p1=poly(X(1);p2=poly(X(2);p3=poly(X(3);p4=poly(X(4);第30页/共160页l01=conv(conv(p2,p3),p4)/(X(1)-X(2)*(X(1)-X(3)*(X(1)-X(4),l11=conv(conv(p1,p3),p4)/(X(2)-X(1)*(X(2)-X(3)*(X(2)-X(4),l21=conv(conv(p1,p2),p4)/(X(3)-X(1)*(X(3)-X(2)*(X(3)-X(4),l31=conv(conv(p1,p2),p3)/(X(4)-X(1)*(X(4)-X(2)*(X(4)-X(3),l0=poly2sym(l01),l1=poly2sym(l11),l2=poly2sym(l21),l3=poly2sym(l31),P=l01*Y(1)+l11*Y(2)+l21*Y(3)+l31*Y(4),运行后输出基函数l0，l1，l2和l3及其插值多项式的系数向量P为l0=-1/24*x3+1/8*x2-1/12*x，l1=1/4*x3-1/4*x2-x+1l2=-1/3*x3+4/3*x，l3=1/8*x3+1/8*x2-1/4*xP=1 4 -4 1第31页/共160页输入程序 L=poly2sym(P),x=0.6;Y=polyval(P,x)运行后输出插值多项式和插值为L=x3+4*x2-4*x+1 Y=0.2560.输入程序 syms M;x=0.6;R3=M*abs(x-X(1)*(x-X(2)*(x-X(3)*(x-X(4)/24运行后输出误差限为R3=91/2500*M 即syms M为定义一个符号变量M四、拉格朗日多项式和基函数的四、拉格朗日多项式和基函数的Matlab程序程序求拉格朗日插值多项式和基函数的求拉格朗日插值多项式和基函数的MATLAB主程序如下主程序如下第32页/共160页function C,L,L1,l=lagran1(X,Y)m=length(X);L=ones(m,m);for k=1:m V=1;for i=1:m if k=i V=conv(V,poly(X(i)/(X(k)-X(i);endendL1(k,:)=V;l(k,:)=poly2sym(V)endC=Y*L1;L=Y*lreturn新建一个lagran1.m文件，将程序拷贝到文件中，并将文件保存到matlab的工作目录下，即可直接调用函数。第33页/共160页例例*1-5给出节点数据f(-2.15)=17.03，f(-1.00)=7.24，f(0.01)=1.05，f(1.02)=2.03，f(2.03)=17.06，f(3.25)=23.05，作五次拉格朗日插值多项式和基函数，并写出估计其误差的公式。解解：输入程序 X=-2.15 -1.00 0.01 1.02 2.03 3.25;Y=17.03 7.24 1.05 2.03 17.06 23.05;C,L,L1,l=lagran1(X,Y)运行后输出五次拉格朗日插值多项式L及其系数向量C，基函数l及其系数矩阵L1如下C=-0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L=1.0954-4.5745*x+3.3960*x2+2.1076*x3+0.0648*x4-0.2169*x5第34页/共160页L1=-0.0056 0.0299 -0.0323 -0.0292 0.0382 -0.0004 0.0331 -0.1377 -0.0503 0.6305 -0.4852 0.0048 -0.0693 0.2184 0.3961 -1.2116 -0.3166 1.0033 0.0687 -0.1469 -0.5398 0.6528 0.9673 -0.0097 -0.0317 0.0358 0.2530 -0.0426 -0.2257 0.0023 0.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002l=-0.0056*x5+0.0299*x4-0.0323*x3-0.0292*x2+0.0382*x-0.0004 0.0331*x5-0.1377*x4-0.0503*x3+0.6305*x2-0.4852*x+0.0048 -0.0693*x5+0.2184*x4+0.3961*x3-1.2116*x2-0.3166*x+1.0033 0.0687*x5-0.1469*x4-0.5398*x3+0.6528*x2+0.9673*x-0.0097 -0.0317*x5+0.0358*x4+0.2530*x3-0.0426*x2-0.2257*x+0.0023 0.0049*x5+0.0004*x4-0.0266*x3+0.0001*x2+0.0220*x-0.0002估计其误差的公式为第35页/共160页五、拉格朗日插值及其误差估计的五、拉格朗日插值及其误差估计的Matlab程序程序function y,R=lagranzi(X,Y,x,M)n=length(X);m=length(x);for i=1:m z=x(i);s=0.0;for k=1:n p=1.0;q1=1.0;c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j);end q1=abs(q1*(z-X(j);c1=c1*j;end s=p*Y(k)+s;end y(i)=s;endR=M*q1/c1;return第36页/共160页例例*1-6已知sin30 =0.5，sin45 =0.7071，sin60 =0.8660，用拉格朗日及其误差估计的MATLAB函数求sin40 的近似值，并估计其误差。解解：输入程序 x=2*pi/9;M=1;X=pi/6,pi/4,pi/3;Y=0.5,0.7071,0.8660;y,R=lagranzi(X,Y,x,M)运行后输出插值y及其误差限R为y=R=0.6434 8.8610e-004.sin(40)=0.6427876第37页/共160页)()(0010101xxxxyyyxP +=)(001010 xxxxxxy +=()fffx0,x1 二次牛顿插值多项式二次牛顿插值多项式 我们再看线性插值线性插值的点斜式点斜式:)(00 xxy+=fx0,x1常数常数(差商差商)由此启发，我们希望二次插值也能类似地有有规律的由此启发，我们希望二次插值也能类似地有有规律的组合表达式组合表达式:P2(x)=0+1(x-x0)+2(x-x0)(x-x1)利用利用P2(x0)=y0有有:0=y0,利用利用P2(x1)=y1有有:1=0101xxxx ()ff=fx0,x1,利用利用P2(x2)=y2有有:2=fx0,x1(x2-x0)(x2-x1)(x2-x0)(x2-x1)0 xx2()ff(x2-x0)-fx0,x2fx0,x1 x2-x1=-=fx0,x1,x2;P2(x)=f(x0)+(x-x0)+(x-x0)(x-x1)fx0,x1 fx0,x1,x2 fx0,x2 x=x0时0第38页/共160页注注:1.事实上事实上,从上述可看出二次牛顿插值公式是用从上述可看出二次牛顿插值公式是用待定系数法待定系数法求得的求得的;2.它也可看作是三个特殊函数的一种线性组合它也可看作是三个特殊函数的一种线性组合:P2(x)=f(x0)+(x-x0)+(x-x0)(x-x1)fx0,x1 fx0,x1,x2 fx0,x1,fx0,x1,x2 f(x0),1,(x-x0),(x-x0)(x-x1)即函数 的线性组合,组合系数为 本质上还是本质上还是基函数法基函数法.第39页/共160页 更一般地，更一般地，n+1个节点的插值多项式，我们希望由上述类似的一组特殊函数：个节点的插值多项式，我们希望由上述类似的一组特殊函数：来线性组合为：来线性组合为：1,(x-x0),(x-x0)(x-x1)，(x-x0)(x-x1)(x-xn)那么其组合系数是什么样的呢？怎么求呢？那么其组合系数是什么样的呢？怎么求呢？我们同样可用待定系数法我们同样可用待定系数法.容易发现容易发现,计算计算a0,a1,a2,an 是很有规律的是很有规律的.第40页/共160页一、均差及其性质一、均差及其性质2 牛顿插值牛顿插值当当x=x0时，时，Pn(x0)=a0=f0.当当x=x1时，时，Pn(x1)=a0+a1(x1-x0)=f1,推推得得a1=f1-f0 x1-x0当当x=x2时，时，Pn(x2)=a0+a1(x2-x0)+a2(x2-x0)(x2-x1)=f2,推得推得f2-f0 x2-x0-f1-f0 x1-x0a2=x2-x1 依次递推可得到依次递推可得到a3,an.为写出系数为写出系数 ak的一般表达式的一般表达式,先引进如下均差定义先引进如下均差定义.第41页/共160页 定义定义2 称称 为函数为函数f(x)关于点关于点x0,xk的的一阶一阶均差均差.称称 为为f(x)的的二阶均差二阶均差.一般地一般地,称称 为为 f(x)的的k 阶均差阶均差(差商差商).fx0,xk=f(xk)-f(x0)xk-x0 fx0,x1,xk=fx0,xk-fx0,x1xk-x1 fx0,x1,xk=fx0,xk-2,xk-fx0,x1,xk-1xk-xk-1均差有如下的基本性质均差有如下的基本性质:1 k 阶均差可表示为函数值阶均差可表示为函数值f(x0),f(x1),f(xk)的线性组合的线性组合,即即 fx0,x1,xk=f(xj)(xj-xj+1)(xj-xk)(xj-xj-1)(xj-x0)kj=0这个性质可用归纳法证明这个性质可用归纳法证明.这个性质也表明均差与节点的排列这个性质也表明均差与节点的排列次序无关次序无关,称为均差的对称性称为均差的对称性,即即 fx0,x1,xk=fx1,x0,x2,xk=fx1,xk,x0第42页/共160页 fx0,x1,xk=f(n)()n!3 若若f(x)在在a,b上存在上存在n阶导数阶导数,且节点且节点x0,x1,xn a,b,则则n阶均差与导数关系如下阶均差与导数关系如下:这个公式可直接用罗尔定理证明这个公式可直接用罗尔定理证明.所以所以 fx0,x1,xk=fx1,xk-1,xk-fx0,x1,xk-1xk-x02 由性质由性质1可得可得:第43页/共160页12 n 1依次将后一式代入前一式依次将后一式代入前一式Nn(x)Rn(x)ai=f x0,xi 二、牛顿插值公式二、牛顿插值公式第44页/共160页Rn(x)Nn(x)n+1(x)多项式多项式Nn(x)显然满足插值条件显然满足插值条件,即即Nn(xj)=f(xj),(j=1,n),且次数不超过且次数不超过n,由唯一性定理它就是前述的由唯一性定理它就是前述的Ln(x),其系数为其系数为 Nn(x)称为牛顿均差插值多项式称为牛顿均差插值多项式,它比拉格朗日插值多项式它比拉格朗日插值多项式计算量省计算量省,且便于程序设计且便于程序设计.第45页/共160页注：注：由由唯一性可知唯一性可知 Nn(x)Ln(x)，只是算法不同，故其余项也相同，即只是算法不同，故其余项也相同，即 实际计算过程为实际计算过程为f(x0)f(x1)f(x2)f(xn 1)f(xn)f x0,x1f x1,x2 f xn 1,xnf x0,x1,x2 f xn 2,xn 1,xnf x0,xn f(xn+1)f xn,xn+1 f xn 1,xn,xn+1 f x1,xn+1 f x0,xn+1均差计算可列均差表如下：均差计算可列均差表如下：第46页/共160页 例例1 依据如下函数值表建立不超过依据如下函数值表建立不超过3次的拉格朗日插值多次的拉格朗日插值多项式及牛顿插值多项式项式及牛顿插值多项式Nn(x),并验证插值多项式的唯一性并验证插值多项式的唯一性.解解:(1)拉格朗日插值多项式拉格朗日插值多项式Ln(x).插值基函数插值基函数xk0124f(xk)19233拉格朗日插值多项式为拉格朗日插值多项式为:第47页/共160页xkf(xk)一阶均差一阶均差二阶均差二阶均差三阶均差三阶均差0119822314343-10-8(2)牛顿插值多项式牛顿插值多项式Nn(x).建立如下差商表建立如下差商表牛顿插值多项式为牛顿插值多项式为:(3)唯一性验证唯一性验证.通过比较牛顿插值多项式和拉格朗日插值多项式通过比较牛顿插值多项式和拉格朗日插值多项式,知知:Nn(x)=Ln(x)这一事实与插值多项式的唯一性一致这一事实与插值多项式的唯一性一致.第48页/共160页例例2注意：当题目中没有指明用哪一种方法建立插值多项式时，原注意：当题目中没有指明用哪一种方法建立插值多项式时，原则上拉格朗日和牛顿插值法都可行。近似计算时，牛顿插值法则上拉格朗日和牛顿插值法都可行。近似计算时，牛顿插值法计算量小，当增加节点时只需增加一项，因而比较方便。相对计算量小，当增加节点时只需增加一项，因而比较方便。相对之下，拉格朗日插值法没有上述优点，但它在理论证明上因其之下，拉格朗日插值法没有上述优点，但它在理论证明上因其插值基函数的许多特点而得到广泛应用。插值基函数的许多特点而得到广泛应用。已知函数已知函数y=f(x)的数据如下表的数据如下表i0123xi0123yi=f(xi)13927试作一个三次插值多项式试作一个三次插值多项式P3(x)计算计算第49页/共160页解：解：利用牛顿插值公式可得利用牛顿插值公式可得为此先作差商表如下为此先作差商表如下kxkf(xk)差商差商00111322296233271864/3第50页/共160页故故由差商表知，由差商表知，而而故令故令即得即得第51页/共160页 等距节点插值公式等距节点插值公式 向前差分向前差分 iiifff=+1ikikikikffff1111)(+=向后差分向后差分 111 =ikikikfffi 1iifff=中心差分中心差分 其中其中当节点当节点等距等距分布时分布时:第52页/共160页 差分的重要性质：差分的重要性质：线性：例如线性：例如 若若 f(x)是是 m 次多项式，则次多项式，则 是是 次多项式，次多项式，而而 差分值可由函数值算出：差分值可由函数值算出：=+=njjknjknfjnf0)1(=+=njnjkjnknfjnf0)1(其中其中 函数值可由差分值算出：函数值可由差分值算出：kjnjknfjnf =+=0kkkhkfxxf!,.,00=knkknnnhkfxxxf!,.,1=kkkhff0)()(=由由 Rn 表达式表达式第53页/共160页牛顿公式牛顿公式 牛顿前差公式牛顿前差公式 牛顿后差公式牛顿后差公式 将节点顺序倒置：将节点顺序倒置：设设，则，则)()()(000 xfkthtxNxNknknn=+=设设，则，则)()1()()(0nknkknnnxfkthtxNxN =+=注：注：一般当一般当 x 靠近靠近 x0 时用前插，靠近时用前插，靠近 xn 时用后插，故两种公式亦称为时用后插，故两种公式亦称为表初公式表初公式和和表末公式表末公式。第54页/共160页*2 牛顿插值的牛顿插值的Matlab实现实现 拉格朗日插值的优点是格式整齐和规范，有误拉格朗日插值的优点是格式整齐和规范，有误差估计公式，它的缺点是没有承袭性，当需要增加差估计公式，它的缺点是没有承袭性，当需要增加节点时，必须重新计算插值的基函数。本节给出具节点时，必须重新计算插值的基函数。本节给出具有承袭性的牛顿插值法的有承袭性的牛顿插值法的MATLAB实现程序。作实现程序。作为基础，首先介绍与牛顿插值有关的差商等计算的为基础，首先介绍与牛顿插值有关的差商等计算的实现程序。实现程序。第55页/共160页一、牛顿插值多项式、差商和误差公式的一、牛顿插值多项式、差商和误差公式的Matlab程序程序输入参数：输入参数：n+1个节点个节点(xi,yi)(i=1,2,n+1)的横坐标向量的横坐标向量X和纵坐标向量和纵坐标向量Y输出参数：输出参数：n阶牛顿插值多项式阶牛顿插值多项式L及其系数向量及其系数向量C，差商的矩阵，差商的矩阵A，插值余项，插值余项公式公式wcgs及其及其 的系数向量的系数向量CW将程序保存为将程序保存为newpoly.m的的M文件，保存到文件，保存到MATLAB工作目录工作目录下即可使用函数。下即