《数学模型与数学实验》课程书籍文献（数学建模算法大全）第09章插值与拟合.pdf_大学文库

-175- 第九章插值与拟合插值：求过已知有限个数据点的近似函数。拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小。插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。 §1 插值方法下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite 插值和三次样条插值。 1.1 拉格朗日多项式插值 1.1.1 插值多项式用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数 f (x)在区间[a,b]上n +1个不同点 n x , x , , x 0 1 L 处的函数值 ( ) i i y = f x (i = 0,1,L,n) ，求一个至多n 次多项式 n n n x = a + a x +L+ a x 0 1 ϕ ( ) （1）使其在给定点处与 f (x)同值，即满足插值条件 (x ) f (x ) y (i 0,1, ,n) ϕ n i = i = i = L （2） (x) ϕ n 称为插值多项式，x (i 0,1, ,n) i = L 称为插值节点，简称节点，[a,b]称为插值区间。从几何上看，n 次多项式插值就是过n +1个点( , ( )) i i x f x (i = 0,1,L,n) ，作一条多项式曲线 y (x) = ϕ n 近似曲线 y = f (x) 。 n 次多项式（1）有n +1个待定系数，由插值条件（2）恰好给出n +1个方程 ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + + + + = + + + + = + + + + = n n n n n n n n n n a a x a x a x y a a x a x a x y a a x a x a x y L LLLLLLLLLLLL L L 2 0 1 2 1 1 2 0 1 1 2 1 0 0 2 0 1 0 2 0 （3）记此方程组的系数矩阵为 A ，则 n n n n n n x x x x x x x x x A L LLLLLLL L L 2 1 2 1 1 0 2 0 0 1 1 1 det( ) = 是范德蒙特(Vandermonde)行列式。当 n x , x , , x 0 1 L 互不相同时，此行列式值不为零。因此方程组（3）有唯一解。这表明，只要n +1个节点互不相同，满足插值要求（2）的插值多项式（1）是唯一的。插值多项式与被插函数之间的差 R (x) f (x) (x) n = −ϕ n

-178- ( ) [ , , , , ] ( ) ( )( ) ( ) [ , , , , ] ( )( ) ( ) [ , , , ] ( ) ( ) ( ) [ , ] 1 0 1 0 1 0 1 0 1 1 0 1 0 0 0 1 n n n n n n n n x f x x x x R x x x x x x x f x x x x x x x x x x f x x x N x f x x x f x x L L L L L L + − = = − − − + − − − = + − + ω 显然，N (x) n 是至多n 次的多项式，且满足插值条件，因而它是 f (x)的n 次插值多项式。这种形式的插值多项式称为 Newton 插值多项式。 R (x) n 称为 Newton 插值余项。 Newton 插值的优点是：每增加一个节点，插值多项式只增加一项，即 ( ) ( ) ( ) ( ) [ , , , ] n+1 = n + − 0 − n 0 1 n+1 N x N x x x L x x f x x L x 因而便于递推运算。而且 Newton 插值的计算量小于 Lagrange 插值。由插值多项式的唯一性可知，Newton 插值余项与 Lagrange 余项也是相等的，即 ( ) ( , ) ( 1)! ( ) ( ) ( ) [ , , , , ] 1 ( 1) 1 0 1 x a b n f R x x f x x x x n n n n n ∈ + = = + + + ω ξ ξ ω L 由此可得差商与导数的关系 ! ( ) [ , , , ] ( ) 0 1 n f f x x x n n ξ L = 其中 ( , ), min{ }, max{ } 0 0 i i n i i n x x ≤ ≤ ≤ ≤ ξ ∈ α β α = β = 。 1.2.3 差分当节点等距时，即相邻两个节点之差（称为步长）为常数，Newton 插值公式的形式会更简单。此时关于节点间函数的平均变化率（差商）可用函数值之差（差分）来表示。定义设有等距节点 ( 0,1, , ) xk = x0 + kh k = L n ，步长 h 为常数， ( ) k k f = f x 。称相邻两个节点 1 , k k + x x 处的函数值的增量 ( 0,1, , 1) f k +1 − f k k = L n − 为函数 f (x) 在点 k x 处以h 为步长的一阶差分，记为 k Δf ，即 ( 0,1, , ) Δf k = f k +1 − f k k = L n 类似地，定义差分的差分为高阶差分。如二阶差分为 ( 0,1, , 2) 1 2 Δ f k = Δf k + − Δf k k = L n − 一般地，m 阶差分为 ( 2,3, ) 1 1 Δm f k = Δm−1 f k + − Δm− f k k = L ，上面定义的各阶差分又称为向前差分。常用的差分还有两种： ∇ k = k − k −1 f f f 称为 f (x)在 k x 处以h 为步长的向后差分； ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ − − ⎠ ⎞ ⎜ ⎝ ⎛ = + 2 2 h f x h f f x δ k k k 称为 f (x) 在 k x 处以 h 为步长的中心差分。一般地， m 阶向后差分与 m 阶中心差分公式为

-179- 1 1 1 − − − ∇ = ∇ − ∇ k m k m k m f f f 2 1 1 2 1 1 − − + − = − k m k m k m δ f δ f δ f 差分具有以下性质：（i）各阶差分均可表成函数值的线性组合，例如 k m j m j j k m f j m f + − = ∑ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Δ = (−1) 0 k j m j j k m f j m f − = ∑ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∇ = (−1) 0 （ii）各种差分之间可以互化。向后差分与中心差分化成向前差分的公式如下： k m m k m f f ∇ = Δ − 2 m k m k m f f − δ = Δ 1.2.4 等距节点插值公式如果插值节点是等距的，则插值公式可用差分表示。设已知节点 xk = x0 + kh (k = 0,1,2,L,n)，则有 ( )( ) ( ) ! ( ) [ , , , ]( )( ) ( ) ( ) ( ) [ , ]( ) 0 1 1 0 0 0 0 0 1 0 1 1 0 0 1 0 − − − − − Δ − + + Δ = + + + − − − = + − n n n n n n x x x x x x n h f x x h f f f x x x x x x x x x N x f x f x x x x L L L L L 若令 x = x0 + th ，则上式又可变形为 0 0 0 0 ! ( 1) ( 1) ( ) f n t t t n N x th f t f n n Δ − − + + = + Δ + + L L 上式称为 Newton 向前插值公式。 1.3 分段线性插值 1.3.1 插值多项式的振荡用 Lagrange 插值多项式 L (x) n 近似 f (x)(a ≤ x ≤ b) ，虽然随着节点个数的增加， L (x) n 的次数n 变大，多数情况下误差| R (x) | n 会变小。但是n 增大时， L (x) n 的光滑性变坏，有时会出现很大的振荡。理论上，当n → ∞ ，在[a,b]内并不能保证 L (x) n 处处收敛于 f (x)。Runge 给出了一个有名的例子： , [ 5,5] 1 1 ( ) 2 ∈ − + = x x f x 对于较大的| x |，随着n 的增大，L (x) n 振荡越来越大，事实上可以证明，仅当| x |≤ 3.63 时，才有 lim L (x) f (x) n n = →∞ ，而在此区间外， L (x) n 是发散的。高次插值多项式的这些缺陷，促使人们转而寻求简单的低次多项式插值。 1.3.2 分段线性插值简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作 I (x) n ，它满足 n i i I (x ) = y ，且 I (x) n 在每个小区间[ , ] i i+1 x x 上是线性

-181- ∑= = − − + n i i i i i i i H x h x x a y y y 0 ( ) [( )(2 ' ) ] 其中 ∏ ≠ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = n j i j i j j i x x x x h 0 2 ， ∑ ≠ = − = n j i j i j i x x a 0 1 。 1.4.2 用 Matlab 实现 Hermite 插值 Matlab 中没有现成的 Hermite 插值函数，必须编写一个 M 文件实现插值。设n 个节点的数据以数组 x0 （已知点的横坐标）， y0 （函数值）， y1（导数值）输入（注意 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）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间[a,b] 的一个分划 Δ : a = x0 < x1 <L< xn−1 < xn = b 如果函数 s(x) 满足：（i）在每个小区间[ , ]( 0,1, , 1) xi xi−1 i = L n − 上 s(x) 是 k 次多项式；（ii） s(x) 在[a,b]上具有 k −1阶连续导数。则称 s(x) 为关于分划 Δ 的k 次样条函数，其图形称为k 次样条曲线。 n x , x , , x 0 1 L 称为样条节点， 1 2 1 , , , n− x x L x 称为内节点， n x , x 0 称为边界点，这类样条函数的全体记做

-182- S ( , k) P Δ ，称为 k 次样条函数空间。显然，折线是一次样条曲线。若 s(x) S ( , k) ∈ P Δ ，则 s(x) 是关于分划 Δ 的 k 次多项式样条函数。k 次多项式样条函数的一般形式为 ∑ ∑ − = + = = + − 1 0 1 ( ) ! ! ( ) n j k j j k i i i k x x i k x s x α β ，其中 (i 0,1, , k) αi = L 和 ( j = 1,2, , n −1) β j L 均为任意常数，而 ⎪⎩ ⎪ ⎨ ⎧ < − ≥ − + = j j k k j j x x x x x x x x 0, ( ) , ( ) ，（ j = 1,2,L, n −1）在实际中最常用的是 k = 2和 3 的情况，即为二次样条函数和三次样条函数。二次样条函数：对于[a,b]上的分划 Δ : a = x0 < x1 <L< xn = b ，则 ∑ − = = + + + − + ∈ Δ 1 1 2 2 2 2 0 1 ( ) ( ,2) 2! 2! ( ) n j j P j s x x x x x S α β α α ，（5）其中 ⎪⎩ ⎪ ⎨ ⎧ < − ≥ − + = j j j j x x x x x x x x 0, ( ) , ( ) 2 2 ，（ j =1,2,L, n −1）。三次样条函数：对于[a,b]上的分划 Δ : a = x0 < x1 <L< xn = b ，则 ∑ − = = + + + + − + ∈ Δ 1 1 2 2 3 3 3 3 0 1 ( ) ( ,3) 2! 3! 3! ( ) n j j P j s x x x x x x S α α β α α ，（6）其中 ⎪⎩ ⎪ ⎨ ⎧ < − ≥ − + = j j j j x x x x x x x x 0, ( ) , ( ) 3 3 ，（ j = 1,2,L, n −1）。利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。例如分段线性插值是一次样条插值。下面我们介绍二次、三次样条插值。 1.5.2 二次样条函数插值首先，我们注意到 ( ) ( ,2) s2 x ∈ SP Δ 中含有 n + 2个特定常数，故应需要 n + 2个插值条件，因此，二次样条插值问题可分为两类：问题（1）：已知插值节点 i x 和相应的函数值 y (i 0,1, , n) i = L 以及端点 0 x （或 n x ）处的导数值 0 y' （或 n y' ），求 ( ) ( ,2) s2 x ∈ SP Δ 使得 ⎩ ⎨ ⎧ = = = = ' ( ) ( ' ( ) ) ( ) ( 0,1,2, , ) 2 0 0 2 n n n i i s x y s x y s x y i n 或 L （7）问题（2）：已知插值节点 i x 和相应的导数值 y' (i 0,1,2, ,n) i = L 以及端点 0 x （或 n x ）处的函数值 0 y （或 n y ），求 ( ) ( ,2) s2 x ∈ SP Δ 使得 ⎩ ⎨ ⎧ = = = = ( ) ( ( ) ) ' ( ) ' ( 0,1,2, , ) 2 0 0 2 2 n n i i s x y s x y s x y i n 或 L （8）

-183- 事实上，可以证明这两类插值问题都是唯一可解的。对于问题（1），由条件（7） ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ = + = = + + + − = = = + + = = + + = ∑ − = 2 0 1 2 0 0 1 1 2 2 2 0 1 2 1 2 2 1 0 1 1 2 1 0 2 2 0 0 1 0 2 0 ' ( ) ' ( ) ( 2,3, , ) 2 1 2 1 ( ) 2 1 ( ) 2 1 ( ) s x x y s x x x x x y j n s x x x y s x x x y j i j j j i j i j α α α α α β α α α α α α L 引入记号 T X n ( , , , , , ) = α 0 α1 α 2 β1 L β −1 为未知向量， ( , , , , ' ) 0 1 0 C y y y y = L n 为已知向量。 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − = − − 0 1 0 0 ( ) 2 1 ( ) 2 1 2 1 1 ( ) 0 2 1 2 1 1 0 0 2 1 1 0 0 2 1 1 0 2 1 2 1 2 2 2 1 2 2 2 2 1 1 2 0 0 L L M M M M M L L L x x x x x x x x x x x x x x x A n n n n n 于是，问题转化为求方程组 AX = C 的解 T X n ( , , , , , ) = α 0 α1 α 2 β1 L β −1 的问题，即可得到二次样条函数 ( ) 2 s x 的表达式。对于问题（2）的情况类似。 1.5.3 三次样条函数插值由于 ( ) ( ,3) s3 x ∈ SP Δ 中含有 n + 3 个待定系数，故应需要 n + 3 个插值条件，已知插值节点 i x 和相应的函数值 f (x ) y (i 0,1,2, , n) i = i = L ，这里提供了 n +1个条件，还需要 2 个边界条件。常用的三次样条函数的边界条件有 3 种类型：（i） n s' (a) y' ,s' (b) y' 3 = 0 3 = 。由这种边界条件建立的样条插值函数称为 f (x)的完备三次样条插值函数。特别地， y' 0 = y' n = 0 时，样条曲线在端点处呈水平状态。如果 f '(x) 不知道，我们可以要求 ' ( ) 3 s x 与 f '(x) 在端点处近似相等。这时以 0 1 2 3 x , x , x , x 为节点作一个三次 Newton 插值多项式 N (x) a ，以 1 2 3 , , , n n− n− n− x x x x 作一个三次 Newton 插值多项式 N (x) b ，要求 s'(a) N' (a), s'(b) N' (b) = a = b 由这种边界条件建立的三次样条称为 f (x)的 Lagrange 三次样条插值函数。（ii） 3 0 3 3 s" (a) = y" ,s" (b) = y" 。特别地 y"n = y"n = 0 时，称为自然边界条件

-184- （iii） ' ( 0) ' ( 0), " ( 0) " ( 0) s 3 a + = s 3 b − s 3 a + = s 3 b − ，(这里要求 s3 (a + 0) = ( 0) s3 b − )此条件称为周期条件。 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 的一个1×2 矩阵来表示，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 机床加工待加工零件的外形根据工艺要求由一组数据(x, y) 给出（在平面情况下），用程控铣床加工时每一刀只能沿 x 方向和 y 方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的(x, y) 坐标。表 1 中给出的 x, y 数据位于机翼断面的下轮廓线上，假设需要得到 x 坐标每改变 0.1 时的 y 坐标。试完成加工所需数据，画出曲线，并求出 x = 0 处的曲线斜率和 13 ≤ x ≤ 15 范围内 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];

《数学模型与数学实验》课程书籍文献（数学建模算法大全）第09章 插值与拟合

《数学模型与数学实验》课程书籍文献（数学建模算法大全）第09章插值与拟合