第09章 插值与拟合

-175- 第九章 插值与拟合 插值：求过已知有限个数据点的近似函数。 拟合：已知有限个数据点，求近似函数，不 要求 对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗 过已知数据点，只要求在某种意义 下它在这些点上的总偏差最小。 插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二 者的数学方法上是完全不同的。而面对一个实际问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 ，究竟应该用插值还是拟合，有时 容易确定，有时则并不明显。 §1 插值方法 下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插 值、Hermite 插值和三次样条插值。 1.1 拉格朗日多项式插值 1.1.1 插值多项式 用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数 )(xf 在 区间 ],[ ba 上 1+n 个不同点 nxxx ,,, 10 L 处的函数值 )( ii xfy = ),,1,0( ni L= ，求一个 至多n次多项式 n nn xaxaax +++= L10)(ϕ （1） 使其在给定点处与 )(xf 同值，即满足插值条件 ),,1,0()()( niyxfx iiin L===ϕ （2） )(xnϕ 称为插值多项式， ),,1,0( nixi L= 称为插值节点，简称节点， ],[ ba 称为插值区 间。从几何上看，n次多项式插值就是过 1+n 个点 ))(,( ii xfx ),,1,0( ni L= ，作一条 多项式曲线 )(xy nϕ= 近似曲线 )(xfy = 。 n次多项式（1）有 1+n 个待定系数，由插值条件（2）恰好给出 1+n 个方程 ⎪⎪⎩ ⎪⎪⎨ ⎧ =++++ =++++ =++++ n n nnnn n n n n yxaxaxaa yxaxaxaa yxaxaxaa L LLLLLLLLLLLL L L 2 210 11 2 12110 00 2 02010 （3） 记此方程组的系数矩阵为 A，则 n nnn n n xxx xxx xxx A L LLLLLLL L L 2 1 2 11 0 2 00 1 1 1 )det( = 是范德蒙特(Vandermonde)行列式。当 nxxx ,,, 10 L 互不相同时，此行列式值不为零。因 此方程组（3）有唯一解。这表明，只要 1+n 个节点互不相同，满足插值要求（2）的 插值多项式（1）是唯一的。 插值多项式与被插函数之间的差 )()()( xxfxR nn ϕ−= ZHENDDONG 高亮 ZHENDDONG 矩形 ZHENDDONG 矩形 ZHENDDONG 矩形 ZHENDDONG 矩形 ZHENDDONG 矩形 ZHENDDONG 高亮 -176- 称为截断误差，又称为插值余项。当 )(xf 充分光滑时， ),(),( )!1( )()()()( 1 )1( bax n fxLxfxR n n nn ∈+=−= + + ξωξ 其中 ∏ = + −= n j jn xxx 0 1 )()(ω 。 1.1.2 拉格朗日插值多项式 实际上比较方便的作法不是解方程（3）求待定系数，而是先构造一组基函数 )())(()( )())(()()( 110 110 niiiiii nii i xxxxxxxx xxxxxxxxxl −−−− −−−−= +− +− LL LL )10( , 0 ,n,,i xx xxn ij j ji j L=− −=∏ ≠= )(xli 是n次多项式，满足 ⎩⎨ ⎧ = ≠= ij ij xl ji 1 0 )( 令 ∑ ∑ ∏ = = ≠= ⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ − −== n i n i n ij j ji j iiin xx xx yxlyxL 0 0 0 )()( （4） 上式称为n次 Lagrange 插值多项式，由方程（3）解的唯一性， 1+n 个节点的n次 Lagrange 插值多项式存在唯一。 1.1.3 用 Matlab 作 Lagrange 插值 Matlab中没有现成的Lagrange插值函数，必须编写一个M文件实现Lagrange插值。 设n个节点数据以数组 0,0 yx 输入（注意 Matlat 的数组下标从 1 开始），m个插值 点以数组 x输入，输出数组 y为m个插值。编写一个名为 lagrange.m 的 M 文件： function y=lagrange(x0,y0,x); n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end -177- 1.2 牛顿（Newton）插值 在导出 Newton 公式前，先介绍公式表示中所需要用到的差商、差分的概念及性质。 1.2.1 差商 定义 设有函数 L,,,),( 210 xxxxf 为一系列互不相等的点，称 ji ji xx xfxf − − )()( )( ji ≠ 为 )(xf 关于点 ji xx , 一阶差商（也称均差）记为 ],[ ji xxf ，即 ji ji ji xx xfxf xxf − −= )()(],[ 称一阶差商的差商 ki kjji xx xxfxxf − − ],[],[ 为 )(xf 关于点 kji xxx ,, 的二阶差商，记为 ],,[ kji xxxf 。一般地，称 k kk xx xxxfxxxf − −− 0 21110 ],,,[],,,[ LL 为 )(xf 关于点 kxxx ,,, 10 L 的 k 阶差商，记为 k kk k xx xxxfxxxfxxxf − −= − 0 21110 10 ],,,[],,,[],,,[ LLL 容易证明，差商具有下述性质： ],[],[ ijji xxfxxf = ],,[],,[],,[ kijjkikji xxxfxxxfxxxf == 1.2.2 Newton 插值公式 线性插值公式可表成 ],[)()()( 10001 xxfxxxfx −+=ϕ 称为一次 Newton 插值多项式。一般地，由各阶差商的定义，依次可得 ],[)()()( 000 xxfxxxfxf −+= ],,[)(],[],[ 101100 xxxfxxxxfxxf −+= ],,,[)(],,[],,[ 210221010 xxxxfxxxxxfxxxf −+= LL ],,,[)(],,,[],,,[ 01010 nnnn xxxfxxxxxfxxxf LLL −+=− 将以上各式分别乘以 ,),)((),(,1 100 Lxxxxxx −−− )())(( 110 −−−− nxxxxxx L ，然 后相加并消去两边相等的部分，即得 ],,,,[)())(( ],,,[)())(( ],[)()()( 1010 10110 1000 nn nn xxxxfxxxxxx xxxfxxxxxx xxfxxxfxf LL LL L −−−+ −−−+ +−+= − 记 -178- ],,,,[)( ],,,,[)())(()( ],,,[)())(( ],[)()()( 101 1010 10110 1000 nn nnn nn n xxxxfx xxxxfxxxxxxxR xxxfxxxxxx xxfxxxfxN L LL LL L + − = −−−= −−−+ +−+= ω 显然， )(xNn 是至多n次的多项式，且满足插值条件，因而它是 )(xf 的n次插值 多项式。这种形式的插值多项式称为 Newton 插值多项式。 )(xRn 称为 Newton 插值余 项。 Newton 插值的优点是：每增加一个节点，插值多项式只增加一项，即 ],,,[)()()()( 11001 ++ −−+= nnnn xxxfxxxxxNxN LL 因而便于递推运算。而且 Newton 插值的计算量小于 Lagrange 插值。 由插值多项式的唯一性可知，Newton 插值余项与 Lagrange 余项也是相等的，即 ),()( )!1( )( ],,,,[)()( 1 )1( 101 bax n f xxxxfxxR n n nnn ∈+= = + + + ξωξ ω L 由此可得差商与导数的关系 ! )(],,,[ )( 10 n fxxxf n n ξ=L 其中 }{max},{min),,( 00 iniini xx ≤≤≤≤ ==∈ βαβαξ 。 1.2.3 差分 当节点等距时，即相邻两个节点之差（称为步长）为常数，Newton 插值公式的形 式会更简单。此时关于节点间函数的平均变化率（差商）可用函数值之差（差分）来表 示。 定义 设有等距节点 ),,1,0(0 nkkhxxk L=+= ，步长 h为常数， )( kk xff = 。 称相邻两个节点 1, +kk xx 处的函数值的增量 )1,,1,0(1 −=−+ nkff kk L 为函数 )(xf 在 点 kx 处以h为步长的一阶差分，记为 kfΔ ，即 ),,1,0(1 nkfff kkk L=−=Δ + 类似地，定义差分的差分为高阶差分。如二阶差分为 )2,,1,0(1 2 −=Δ−Δ=Δ + nkfff kkk L 一般地，m阶差分为 ),3,2(11 1 L=Δ−Δ=Δ −+− kfff kmkmkm ， 上面定义的各阶差分又称为向前差分。常用的差分还有两种： 1−−=∇ kkk fff 称为 )(xf 在 kx 处以h为步长的向后差分； ⎟⎠ ⎞⎜⎝ ⎛ −−⎟⎠ ⎞⎜⎝ ⎛ += 22 hxfhxff kkkδ 称为 )(xf 在 kx 处以 h为步长的中心差分。一般地，m阶向后差分与m阶中心差分公 式为 -179- 1 11 − −− ∇−∇=∇ kmkmkm fff 2 1 1 2 1 1 − − + − −= k m k m k m fff δδδ 差分具有以下性质： （i）各阶差分均可表成函数值的线性组合，例如 jmk m j j k m f j m f −+ = ∑ ⎟⎟⎠ ⎞⎜⎜⎝ ⎛−=Δ )1( 0 jk m j j k m f j m f − = ∑ ⎟⎟⎠ ⎞⎜⎜⎝ ⎛−=∇ )1( 0 （ii）各种差分之间可以互化。向后差分与中心差分化成向前差分的公式如下： mk m k m ff −Δ=∇ 2 mk m k m ff − Δ=δ 1.2.4 等距节点插值公式 如果插值节点是等距的，则插值公式可用差分表示。设已知节点 khxxk += 0 ),,2,1,0( nk L= ，则有 )())(( ! )( )())(](,,,[ )](,[)()( 110 0 0 0 0 11010 0100 − − −−−Δ++−Δ+= −−−++ −+= nn n nn n xxxxxx hn fxx h ff xxxxxxxxxf xxxxfxfxN LL LLL 若令 thxx += 0 ，则上式又可变形为 0000 ! )1()1()( f n ntttftfthxN nn Δ+−−++Δ+=+ LL 上式称为 Newton 向前插值公式。 1.3 分段线性插值 1.3.1 插值多项式的振荡 用 Lagrange 插值多项式 )(xLn 近似 ))(( bxaxf ≤≤ ，虽然随着节点个数的增加， )(xLn 的次数n变大，多数情况下误差 |)(| xRn 会变小。但是n增大时， )(xLn 的光滑 性变坏，有时会出现很大的振荡。理论上，当 ∞→n ，在 ],[ ba 内并不能保证 )(xLn 处 处收敛于 )(xf 。Runge 给出了一个有名的例子： ]5,5[, 1 1)( 2 −∈+= xxxf 对于较大的 || x ，随着n的增大， )(xLn 振荡越来越大，事实上可以证明，仅当 63.3|| ≤x 时，才有 )()(lim xfxLnn =∞→ ，而在此区间外， )(xLn 是发散的。 高次插值多项式的这些缺陷，促使人们转而寻求简单的低次多项式插值。 1.3.2 分段线性插值 简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性 插值函数，记作 )(xIn ，它满足 iin yxI =)( ，且 )(xIn 在每个小区间 ],[ 1+ii xx 上是线性 -180- 函数 ),,1,0( ni L= 。 )(xIn 可以表示为 ∑ = = n i iin xlyxI 0 )()( ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ =∈− − =∈− − = + + + − − − 其它 时舍去 时舍去 ， ）（， ） 0 ],[ 0( ],[ , )( 1 1 1 1 1 1 nixxx xx xx ixxx xx xx xl ii ii i ii ii i i )(xIn 有良好的收敛性，即对于 ],[ bax∈ 有， )()(lim xfxInn =∞→ 。 用 )(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.1 Hermite 插值多项式 如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一 阶、二阶甚至更高阶的导数值，这就是 Hermite 插值问题。本节主要讨论在节点处插值 函数与函数的值及一阶导数值均相等的 Hermite 插值。 设已知函数 )(xfy = 在 1+n 个互异节点 nxxx ,,, 10 L 上的函数值 )( ii xfy = ),,1,0( ni L= 和导数值 )('' ii xfy = ),,1,0( ni L= ，要求一个至多 12 +n 次的多项式 )(xH ，使得 ii yxH =)( ii yxH ')(' = ),,1,0( ni L= 满足上述条件的多项式 )(xH 称为 Hermite 插值多项式。 Hermite 插值多项式为 -181- ∑ = +−−= n i iiiiii yyyaxxhxH 0 ])'2)([()( 其中 ∏ ≠= ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ − −= n ij j ji j i xx xx h 0 2 ， ∑ ≠= −= n ij j ji i xx a 0 1 。 1.4.2 用 Matlab 实现 Hermite 插值 Matlab 中没有现成的 Hermite 插值函数，必须编写一个 M 文件实现插值。 设n个节点的数据以数组 0x （已知点的横坐标）， 0y （函数值）， 1y （导数值） 输入（注意 Matlat 的数组下标从 1 开始），m个插值点以数组 x输入，输出数组 y为m 个插值。编写一个名为 hermite.m 的 M 文件： function y=hermite(x0,y0,y1,x); n=length(x0);m=length(x); for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j~=i h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; a=1/(x0(i)-x0(j))+a; end end yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); end y(k)=yy; end 1.5 样条插值 许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外 形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续， 而且要有连续的曲率，这就导致了样条插值的产生。 1.5.1 样条函数的概念 所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木 条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并 使连接点处有连续的曲率。 数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间 ],[ ba 的一个分划 bxxxxa nn =<<<<=Δ −110: L 如果函数 )(xs 满足： （i）在每个小区间 )1,,1,0](,[ 1 −=− nixx ii L 上 )(xs 是 k 次多项式； （ii） )(xs 在 ],[ ba 上具有 1−k 阶连续导数。 则称 )(xs 为关于分划Δ的 k 次样条函数，其图形称为 k 次样条曲线。 nxxx ,,, 10 L 称为 样条节点， 121 ,,, −nxxx L 称为内节点， nxx ,0 称为边界点，这类样条函数的全体记做 -182- ),( kSP Δ ，称为 k 次样条函数空间。 显然，折线是一次样条曲线。 若 ),()( kSxs P Δ∈ ，则 )(xs 是关于分划Δ的 k 次多项式样条函数。 k 次多项式样 条函数的一般形式为 ∑∑ − = + = −+= 1 10 )( !! )( n j k j j k i i i k xxki xxs βα ， 其中 ),,1,0( kii L=α 和 )1,,2,1( −= njj Lβ 均为任意常数，而 ⎪⎩ ⎪⎨ ⎧ < ≥−=− + j j k jk j xx xxxx xx ,0 ,)( )( ，（ 1,,2,1 −= nj L ） 在实际中最常用的是 2=k 和 3 的情况，即为二次样条函数和三次样条函数。 二次样条函数：对于 ],[ ba 上的分划 bxxxa n =<<<=Δ L10: ，则 ∑− = + Δ∈−+++= 1 1 222 102 )2,()(!2!2 )( n j Pj j Sxxxxxs βααα ， （5） 其中 ⎪⎩ ⎪⎨ ⎧ < ≥−=− + j jj j xx xxxx xx ,0 ,)( )( 2 2 ，（ 1,,2,1 −= nj L ）。 三次样条函数：对于 ],[ ba 上的分划 bxxxa n =<<<=Δ L10: ，则 ∑− = + Δ∈−++++= 1 1 33322 103 )3,()(!3!3!2 )( n j Pj j Sxxxxxxs βαααα ， （6） 其中 ⎪⎩ ⎪⎨ ⎧ < ≥−=− + j jj j xx xxxx xx ,0 ,)( )( 3 3 ，（ 1,,2,1 −= nj L ）。 利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。例如分段线性插值 是一次样条插值。下面我们介绍二次、三次样条插值。 1.5.2 二次样条函数插值 首先，我们注意到 )2,()(2 Δ∈ PSxs 中含有 2+n 个特定常数，故应需要 2+n 个插 值条件，因此，二次样条插值问题可分为两类： 问题（1）： 已知插值节点 ix 和相应的函数值 ),,1,0( niyi L= 以及端点 0x （或 nx ）处的导数 值 0'y （或 ny' ），求 )2,()(2 Δ∈ PSxs 使得 ⎩⎨ ⎧ == == ))('()(' ),,2,1,0()( 002 2 nnn ii yxsyxs niyxs 或 L （7） 问题（2）： 已知插值节点 ix 和相应的导数值 ),,2,1,0(' niy i L= 以及端点 0x （或 nx ）处的函 数值 0y （或 ny ），求 )2,()(2 Δ∈ PSxs 使得 ⎩⎨ ⎧ == == ))(()( ),,2,1,0(')(' 2002 2 nn ii yxsyxs niyxs 或 L （8） -183- 事实上，可以证明这两类插值问题都是唯一可解的。 对于问题（1），由条件（7） ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ =+= ==−+++= =++= =++= ∑− = 002102 1 1 22 2102 1 2 1211012 0 2 0201002 ')(' ),,3,2()( 2 1 2 1)( 2 1)( 2 1)( yxxs njyxxxxxs yxxxs yxxxs j i jijijjj αα βααα ααα ααα L 引入记号 TnX ),,,,,( 11210 −= ββααα L 为未知向量， )',,,,( 010 yyyyC nL= 为已 知向量。 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ −− −= − 0010 )( 2 1)( 2 1 2 11 0)( 2 1 2 11 00 2 11 00 2 11 0 2 1 2 1 2 2 12 2 22 2 11 2 00 L L MMMMM L L L x xxxxxx xxxx xx xx A nnnnn 于是，问题转化为求方程组 CAX = 的解 TnX ),,,,,( 11210 −= ββααα L 的问题， 即可得到二次样条函数 )(2 xs 的表达式。 对于问题（2）的情况类似。 1.5.3 三次样条函数插值 由于 )3,()(3 Δ∈ PSxs 中含有 3+n 个待定系数，故应需要 3+n 个插值条件，已知 插值节点 ix 和相应的函数值 ),,2,1,0()( niyxf ii L== ，这里提供了 1+n 个条件，还 需要 2 个边界条件。 常用的三次样条函数的边界条件有 3 种类型： （i） nybsyas ')(',')(' 303 == 。由这种边界条件建立的样条插值函数称为 )(xf 的 完备三次样条插值函数。 特别地， 0''0 == nyy 时，样条曲线在端点处呈水平状态。 如果 )(' xf 不知道，我们可以要求 )('3 xs 与 )(' xf 在端点处近似相等。这时以 3210 ,,, xxxx 为节点作一个三次 Newton 插值多项式 )(xNa ，以 321 ,,, −−− nnnn xxxx 作一 个三次 Newton 插值多项式 )(xNb ，要求 )(')('),(')(' bNbsaNas ba == 由这种边界条件建立的三次样条称为 )(xf 的 Lagrange 三次样条插值函数。 （ii） 3303 ")(",")(" ybsyas == 。特别地 0"" == nn yy 时，称为自然边界条件。 -184- （iii） )0(")0("),0(')0(' 3333 −=+−=+ bsasbsas ，(这里要求 =+ )0(3 as )0(3 −bs )此条件称为周期条件。 1.5.4 三次样条插值在 Matlab 中的实现 在 Matlab 中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法， 就是采用非扭结（not-a-knot）条件。这个条件强迫第 1 个和第 2 个三次多项式的三阶 导数相等。对最后一个和倒数第 2 个三次多项式也做同样地处理。 Matlab 中三次样条插值也有现成的函数： y=interp1(x0,y0,x,'spline')； y=spline(x0,y0,x)； pp=csape(x0,y0,conds)，y=ppval(pp,x)。 其中 x0,y0 是已知数据点，x 是插值点，y 是插值点的函数值。 对于三次样条插值，我们提倡使用函数 csape，csape 的返回值是 pp 形式，要求插 值点的函数值，必须调用函数 ppval。 pp=csape(x0,y0)：使用默认的边界条件，即 Lagrange 边界条件。 pp=csape(x0,y0,conds)中的 conds 指定插值的边界条件，其值可为： 'complete' 边界为一阶导数，即默认的边界条件 'not-a-knot' 非扭结条件 'periodic' 周期条件 'second' 边界为二阶导数，二阶导数的值[0, 0]。 'variational' 设置边界的二阶导数值为[0,0]。 对于一些特殊的边界条件，可以通过 conds 的一个 21× 矩阵来表示，conds 元素的 取值为 1，2。此时，使用命令 pp=csape(x0,y0_ext,conds) 其中 y0_ext=[left, y0, right]，这里 left 表示左边界的取值，right 表示右边界的取值。 conds(i)=j 的含义是给定端点 i的 j阶导数，即 conds 的第一个元素表示左边界的条 件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶 导数，对应的值由 left 和 right 给出。 详细情况请使用帮助 help csape。 例 1 机床加工 待加工零件的外形根据工艺要求由一组数据 ),( yx 给出（在平面情况下），用程控 铣床加工时每一刀只能沿 x方向和 y方向走非常小的一步，这就需要从已知数据得到加 工所要求的步长很小的 ),( yx 坐标。 表 1 中给出的 yx, 数据位于机翼断面的下轮廓线上，假设需要得到 x坐标每改变 0.1 时的 y坐标。试完成加工所需数据，画出曲线，并求出 0=x 处的曲线斜率和 1513 ≤≤ x 范围内 y的最小值。 表 1 x 0 3 5 7 9 11 12 13 14 15 y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 要求用 Lagrange、分段线性和三次样条三种插值方法计算。 解 编写以下程序： clc,clear x0=[0 3 5 7 9 11 12 13 14 15]; -185- y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6]; x=0:0.1:15; y1=lagrange(x0,y0,x); %调用前面编写的Lagrange插值函数 y2=interp1(x0,y0,x); y3=interp1(x0,y0,x,'spline'); pp1=csape(x0,y0); y4=ppval(pp1,x); pp2=csape(x0,y0,'second'); y5=ppval(pp2,x); fprintf('比较一下不同插值方法和边界条件的结果:\n') fprintf('x y1 y2 y3 y4 y5\n') xianshi=[x',y1',y2',y3',y4',y5']; fprintf('%f\t%f\t%f\t%f\t%f\t%f\n',xianshi') subplot(2,2,1), plot(x0,y0,'+',x,y1), title('Lagrange') subplot(2,2,2), plot(x0,y0,'+',x,y2), title('Piecewise linear') subplot(2,2,3), plot(x0,y0,'+',x,y3), title('Spline1') subplot(2,2,4), plot(x0,y0,'+',x,y4), title('Spline2') dyx0=ppval(fnder(pp1),x0(1)) %求x=0处的导数 ytemp=y3(131:151); index=find(ytemp==min(ytemp)); xymin=[x(130+index),ytemp(index)] 计算结果略。 可以看出，拉格朗日插值的结果根本不能应用，分段线性插值的光滑性较差（特别 是在 14=x 附近弯曲处），建议选用三次样条插值的结果。 1.6 B 样条函数插值方法 1.6.1 磨光函数 实际中的许多问题，往往是既要求近似函数（曲线或曲面）有足够的光滑性，又要 求与实际函数有相同的凹凸性，一般插值函数和样条函数都不具有这种性质。如果对于 一个特殊函数进行磨光处理生成磨光函数（多项式），则用磨光函数构造出样条函数作 为插值函数，既有足够的光滑性，而且也具有较好的保凹凸性，因此磨光函数在一维插 值（曲线）和二维插值（曲面）问题中有着广泛的应用。 由积分理论可知，对于可积函数通过积分会提高函数的光滑度，因此，我们可以利 用积分方法对函数进行磨光处理。 定义 若 )(xf 为可积函数，对于 0>h ，则称积分 ∫ +−= 2 2 ,1 )( 1)( hx hxh dttf h xf 为 )(xf 的一次磨光函数，h称为磨光宽度。 同样的，可以定义 )(xf 的 k 次磨光函数为 ∫ +− −= 2 2 ,1, )( 1)( hx hx hkhk dttf h xf （ 1>k ） 事实上，磨光函数 )(, xf hk 比 )(xf 的光滑程度要高，且当磨光宽度h很少时 )(, xf hk 很接近于 )(xf 。 1.6.2 等距 B 样条函数 对于任意的函数 )(xf ，定义其步长为 1 的中心差分算子δ 如下： -186- ) 2 1() 2 1()( −−+= xfxfxfδ 在此取 0)( += xxf ，则 00 0 2 1 2 1 ++ + ⎟⎠ ⎞⎜⎝ ⎛ −−⎟⎠ ⎞⎜⎝ ⎛ += xxxδ 是一个单位方波函数（如图 1），记 00 )( +=Ω xx δ 。并取 1=h ，对 )(0 xΩ 进行一次磨光 得 ∫∫ +− ++ + − ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛ −−⎟⎠ ⎞⎜⎝ ⎛ +=Ω=Ω 2 1 2 1 00 2 1 2 1 01 2 1 2 1)()( x x x x dtttdttx +++− + + + −+−+=−= ∫∫ )1(2)1(1 01 0 xxxdttdtt xxxx 显然 )(1 xΩ 是连续的（如图 1）。 图 1 )(0 xΩ 和 )(1 xΩ 的图形 类似地可得到 k 次磨光函数为 ∑+ = + + ⎟⎠ ⎞⎜⎝ ⎛ −++−=Ω 1 0 1 2 1 ! )1()( k j kj kj k j kx k Cx 实际上，可以证明： )(xkΩ 是分段 k 次多项式，且具有 1−k 阶连续导数，其 k 阶 导数有 2+k 个间断点，记为 2 1+−= kjx j （ 1,,2,1,0 += kj L ）。从而可知 )(xkΩ 是 对应于分划 +∞<<<<<−∞Δ +110: kxxx L 的 k 次多项式样条函数，称之为基本样 条函数，简称为 k 次 B 样条。由于样条节点为 2 1+−= kjx j （ 1,,2,1,0 += kj L ）是 等距的，故 )(xkΩ 又称为 k 次等距 B 样条函数。 对于任意函数 )(xf 的 k 次磨光函数，由归纳法可以得到： dttf h tx h xf khk ∫ ∞+∞− − ⎟⎠⎞⎜⎝⎛ −Ω= )(1)( 1, （ 22 hxthx +≤≤− ） 特别地，当 1)( =xf 时，有 11 1 =⎟⎠ ⎞⎜⎝ ⎛ −Ω∫ ∞+∞− − dth txh k ，从而 1)( =Ω∫ +∞ ∞− dxxk ，且当 1≥k 时有递推关系 -187- ⎥⎦ ⎤⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛ −Ω⎟⎠ ⎞⎜⎝ ⎛ −−−⎟⎠ ⎞⎜⎝ ⎛ +Ω⎟⎠ ⎞⎜⎝ ⎛ ++=Ω −− 2 1 2 1 2 1 2 11)( 11 xx kxkx k x kkk 1.6.3 一维等距 B 样条函数插值 等距 B 样条函数与通常的样条有如下的关系： 定理 设有区间 ],[ ba 的均匀分划 jhxx j +=Δ 0: （ nj ,,2,1,0 L= ）， n abh −= ， 则对任意 k 次样条函数 ),()( kSxs Pk Δ∈ 都可以表示为 B 样条函数族 1 0 2 1 −= −=⎭⎬ ⎫ ⎩⎨ ⎧ ⎟⎠ ⎞⎜⎝ ⎛ +−−−Ω nj kj k kj h xx 的线性组合。 根 据 定 理 ， 如 果 已 知 曲 线 上 一 组 点 ),( jj yx ， 其 中 jhxx j += 0 （ njh ,,2,1,0,0 L=> ），则可以构造出一条样条磨光曲线（即为 B 样条函数族的线 性组合） ∑− −= ⎟⎠ ⎞⎜⎝ ⎛ −−Ω= 1 0)( n kj kjk jh xxcxs 其中 jc （ 1,,1, −+−−= nkkj L ）为待定常数。用它来逼近曲线，既有较好的精度， 又有良好的保凸性。 实际中，最常用的是 3=k 的情况，即一般形式为 ∑+ −= ⎟⎠ ⎞⎜⎝ ⎛ −−Ω= 1 1 0 33 )( n j j jh xxcxs 其中 3+n 个待定系数 jc （ 1,,0,1 +−= nj L ）可以由插值条件确定。 对于插值条件 ⎪⎩ ⎪⎨⎧ == == ),0(')(' ),2,1,0()( 3 3 njyxs njyxs jj jj L 有 ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ =−Ω= ==−Ω= =−Ω= ∑ ∑ ∑ + −= + −= + −= n n j jn i n j ji n j j yjnc h xs niyjicxs yjc h xs ')('1)(' ,,2,1,0,)()( ')('1)(' 3 1 1 3 1 1 33 0 1 1 303 L （9） 注意到 )(3 xΩ 的局部非零性及其函数值： 3 2)0(3 =Ω ， 6 1)1(3 =±Ω ，当 2≥x 时 0)(3 =Ω x ；且由 )2 1() 2 1()(' 223 −Ω−+Ω=Ω xxx 知， 0)0('3 =Ω ， 2 1)1('3 m=±Ω ， 当 2≥x 时 0)('3 =Ω x 。则（9）式的每一个方程只有三个非零系数，具体为 -188- ⎪⎩ ⎪⎨ ⎧ =+− ==++ =+− +− +− − nnn iiii hycc niyccc hycc '2 ,,2,1,0,64 '2 11 11 011 L （10） 由方程组（10）容易求解出 jc （ 1,,0,1 +−= nj L ），即可得到三次样条函数 )(3 xs 表达式。 1.6.4 二维等距 B 样条函数插值 设有空间曲面 ),( yxfz = （未知），如果已知二维等距节点 =),( ji yx ),( 00 τjyihx ++ （ 0, >τh ）上的值为 ijz （ mjni ,,1,0;,,1,0 LL == ），则相应的 B 样条磨光曲面的一般形式为 ⎟⎠ ⎞⎜⎝ ⎛ −−Ω⎟⎠ ⎞⎜⎝ ⎛ −−Ω= ∑∑− −= − −= jyyi h xxcyxs l n ki k m lj ij τ 0 1 0 1 ),( 其中 ijc （ mjni ,,1,0;,,1,0 LL == ）为待定常数， lk, 可以取不同值，常用的也是 2, =lk 和3的情形。这是一种具有良好保凸性的光滑曲面（函数），在工程设计中是常 用的，但只能使用于均匀划分或近似均匀划分的情况。 1.7 二维插值 前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。若 节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的 高程（节点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这 些点的高程（插值）。 1.7.1 插值节点为网格节点 已知 nm× 个节点： ),,( ijji zyx （ njmi ,,2,1;,,2,1 LL == ），且 mxx <