重庆工商大学数学建模算法讲义第09章 插值与拟合

-175-第九章插值与拟合插值：求过已知有限个数据点的近似函数。拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小。插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。§1插值方法下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite插值和三次样条插值。1.1拉格朗日多项式插值1.1.1插值多项式用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数)(xf在区间],[ba上1+n个不同点nxxx,,,10L处的函数值)(iixfy=),,1,0(niL=，求一个至多n次多项式nnnxaxaax+++=L10)(ϕ（1）使其在给定点处与)(xf同值，即满足插值条件),,1,0()()(niyxfxiiinL===ϕ（2）)(xnϕ称为插值多项式，),,1,0(nixiL=称为插值节点，简称节点，],[ba称为插值区间。从几何上看，n次多项式插值就是过1+n个点))(,(iixfx),,1,0(niL=，作一条多项式曲线)(xynϕ=近似曲线)(xfy=。n次多项式（1）有1+n个待定系数，由插值条件（2）恰好给出1+n个方程⎪⎪⎩⎪⎪⎨⎧=++++=++++=++++nnnnnnnnnnyxaxaxaayxaxaxaayxaxaxaaLLLLLLLLLLLLLLL22101121211000202010（3）记此方程组的系数矩阵为A，则nnnnnnxxxxxxxxxALLLLLLLLLL212110200111)det(=是范德蒙特(Vandermonde)行列式。当nxxx,,,10L互不相同时，此行列式值不为零。因此方程组（3）有唯一解。这表明，只要1+n个节点互不相同，满足插值要求（2）的插值多项式（1）是唯一的。插值多项式与被插函数之间的差)()()(xxfxRnnϕ−=-176-称为截断误差，又称为插值余项。当)(xf充分光滑时，),(),()!1()()()()(1)1(baxnfxLxfxRnnnn∈+=−=++ξωξ其中∏=+−=njjnxxx01)()(ω。1.1.2拉格朗日插值多项式实际上比较方便的作法不是解方程（3）求待定系数，而是先构造一组基函数)())(()()())(()()(110110niiiiiiniiixxxxxxxxxxxxxxxxxl−−−−−−−−=+−+−LLLL)10(,0,n,,ixxxxnijjjijL=−−=∏≠=)(xli是n次多项式，满足⎩⎨⎧=≠=ijijxlji10)(令∑∑∏==≠=⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛−−==nininijjjijiiinxxxxyxlyxL000)()(（4）上式称为n次Lagrange插值多项式，由方程（3）解的唯一性，1+n个节点的n次Lagrange插值多项式存在唯一。1.1.3用Matlab作Lagrange插值Matlab中没有现成的Lagrange插值函数，必须编写一个M文件实现Lagrange插值。设n个节点数据以数组0,0yx输入（注意Matlat的数组下标从1开始），m个插值点以数组x输入，输出数组y为m个插值。编写一个名为lagrange.m的M文件：functiony=lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end-177-1.2牛顿（Newton）插值在导出Newton公式前，先介绍公式表示中所需要用到的差商、差分的概念及性质。1.2.1差商定义设有函数L,,,),(210xxxxf为一系列互不相等的点，称jijixxxfxf−−)()()(ji≠为)(xf关于点jixx,一阶差商（也称均差）记为],[jixxf，即jijijixxxfxfxxf−−=)()(],[称一阶差商的差商kikjjixxxxfxxf−−],[],[为)(xf关于点kjixxx,,的二阶差商，记为],,[kjixxxf。一般地，称kkkxxxxxfxxxf−−−021110],,,[],,,[LL为)(xf关于点kxxx,,,10L的k阶差商，记为kkkkxxxxxfxxxfxxxf−−=−02111010],,,[],,,[],,,[LLL容易证明，差商具有下述性质：],[],[ijjixxfxxf=],,[],,[],,[kijjkikjixxxfxxxfxxxf==1.2.2Newton插值公式线性插值公式可表成],[)()()(10001xxfxxxfx−+=ϕ称为一次Newton插值多项式。一般地，由各阶差商的定义，依次可得],[)()()(000xxfxxxfxf−+=],,[)(],[],[101100xxxfxxxxfxxf−+=],,,[)(],,[],,[210221010xxxxfxxxxxfxxxf−+=LL],,,[)(],,,[],,,[01010nnnnxxxfxxxxxfxxxfLLL−+=−将以上各式分别乘以,),)((),(,1100Lxxxxxx−−−)())((110−−−−nxxxxxxL，然后相加并消去两边相等的部分，即得],,,,[)())((],,,[)())((],[)()()(1010101101000nnnnxxxxfxxxxxxxxxfxxxxxxxxfxxxfxfLLLLL−−−+−−−++−+=−记-178-],,,,[)(],,,,[)())(()(],,,[)())((],[)()()(1011010101101000nnnnnnnnxxxxfxxxxxfxxxxxxxRxxxfxxxxxxxxfxxxfxNLLLLLL+−=−−−=−−−++−+=ω显然，)(xNn是至多n次的多项式，且满足插值条件，因而它是)(xf的n次插值多项式。这种形式的插值多项式称为Newton插值多项式。)(xRn称为Newton插值余项。Newton插值的优点是：每增加一个节点，插值多项式只增加一项，即],,,[)()()()(11001++−−+=nnnnxxxfxxxxxNxNLL因而便于递推运算。而且Newton插值的计算量小于Lagrange插值。由插值多项式的唯一性可知，Newton插值余项与Lagrange余项也是相等的，即),()()!1()(],,,,[)()(1)1(101baxnfxxxxfxxRnnnnn∈+==+++ξωξωL由此可得差商与导数的关系!)(],,,[)(10nfxxxfnnξ=L其中}{max},{min),,(00iniinixx≤≤≤≤==∈βαβαξ。1.2.3差分当节点等距时，即相邻两个节点之差（称为步长）为常数，Newton插值公式的形式会更简单。此时关于节点间函数的平均变化率（差商）可用函数值之差（差分）来表示。定义设有等距节点),,1,0(0nkkhxxkL=+=，步长h为常数，)(kkxff=。称相邻两个节点1,+kkxx处的函数值的增量)1,,1,0(1−=−+nkffkkL为函数)(xf在点kx处以h为步长的一阶差分，记为kfΔ，即),,1,0(1nkfffkkkL=−=Δ+类似地，定义差分的差分为高阶差分。如二阶差分为)2,,1,0(12−=Δ−Δ=Δ+nkfffkkkL一般地，m阶差分为),3,2(111L=Δ−Δ=Δ−+−kfffkmkmkm，上面定义的各阶差分又称为向前差分。常用的差分还有两种：1−−=∇kkkfff称为)(xf在kx处以h为步长的向后差分；⎟⎠⎞⎜⎝⎛−−⎟⎠⎞⎜⎝⎛+=22hxfhxffkkkδ称为)(xf在kx处以h为步长的中心差分。一般地，m阶向后差分与m阶中心差分公式为-179-111−−−∇−∇=∇kmkmkmfff211211−−+−−=kmkmkmfffδδδ差分具有以下性质：（i）各阶差分均可表成函数值的线性组合，例如jmkmjjkmfjmf−+=∑⎟⎟⎠⎞⎜⎜⎝⎛−=Δ)1(0jkmjjkmfjmf−=∑⎟⎟⎠⎞⎜⎜⎝⎛−=∇)1(0（ii）各种差分之间可以互化。向后差分与中心差分化成向前差分的公式如下：mkmkmff−Δ=∇2mkmkmff−Δ=δ1.2.4等距节点插值公式如果插值节点是等距的，则插值公式可用差分表示。设已知节点khxxk+=0),,2,1,0(nkL=，则有)())((!)()())(](,,,[)](,[)()(1100000110100100−−−−−Δ++−Δ+=−−−++−+=nnnnnnxxxxxxhnfxxhffxxxxxxxxxfxxxxfxfxNLLLLL若令thxx+=0，则上式又可变形为0000!)1()1()(fnntttftfthxNnnΔ+−−++Δ+=+LL上式称为Newton向前插值公式。1.3分段线性插值1.3.1插值多项式的振荡用Lagrange插值多项式)(xLn近似))((bxaxf≤≤，虽然随着节点个数的增加，)(xLn的次数n变大，多数情况下误差|)(|xRn会变小。但是n增大时，)(xLn的光滑性变坏，有时会出现很大的振荡。理论上，当∞→n，在],[ba内并不能保证)(xLn处处收敛于)(xf。Runge给出了一个有名的例子：]5,5[,11)(2−∈+=xxxf对于较大的||x，随着n的增大，)(xLn振荡越来越大，事实上可以证明，仅当63.3||≤x时，才有)()(limxfxLnn=∞→，而在此区间外，)(xLn是发散的。高次插值多项式的这些缺陷，促使人们转而寻求简单的低次多项式插值。1.3.2分段线性插值简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作)(xIn，它满足iinyxI=)(，且)(xIn在每个小区间],[1+iixx上是线性-180-函数),,1,0(niL=。)(xIn可以表示为∑==niiinxlyxI0)()(⎪⎪⎪⎩⎪⎪⎪⎨⎧=∈−−=∈−−=+++−−−其它时舍去时舍去，）（，）0],[0(],[,)(111111nixxxxxxxixxxxxxxxliiiiiiiiiii)(xIn有良好的收敛性，即对于],[bax∈有，)()(limxfxInn=∞→。用)(xIn计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n无关。但n越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。1.3.3用Matlab实现分段线性插值用Matlab实现分段线性插值不需要编制函数程序，Matlab中有现成的一维插值函数interp1。y=interp1(x0,y0,x,'method')method指定插值的方法，默认为线性插值。其值可为：'nearest'最近项插值'linear'线性插值'spline'逐段3次样条插值'cubic'保凹凸性3次插值。所有的插值方法要求x0是单调的。当x0为等距时可以用快速插值法，使用快速插值法的格式为'*nearest'、'*linear'、'*spline'、'*cubic'。1.4埃尔米特(Hermite)插值1.4.1Hermite插值多项式如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值，这就是Hermite插值问题。本节主要讨论在节点处插值函数与函数的值及一阶导数值均相等的Hermite插值。设已知函数)(xfy=在1+n个互异节点nxxx,,,10L上的函数值)(iixfy=),,1,0(niL=和导数值)(''iixfy=),,1,0(niL=，要求一个至多12+n次的多项式)(xH，使得iiyxH=)(iiyxH'