第九章插值与拟合 插值:求过已知有限个数据点的近似函数。 合:已知有限个数据点,求近似函数,不要求过已知数据点,只要求在某种意义 这些点上的总偏差最小 插值和拟合都是要根据一组数据构造一个函数作为近似,由于近似的要求不同,二 者的数学方法上是完全不同的。而面对一个实际问题,究竟应该用插值还是拟合,有时 容易确定,有时则并不明显。 §1插值方法 下面介绍几种基本的、常用的插值:拉格朗日多项式插值、牛顿插值、分段线性插 值、 Hermite插值和三次样条插值。 1.1拉格朗日多项式插值 1.1.1插值多项式 用多项式作为研究插值的工具,称为代数插值。其基本问题是:已知函数∫(x)在 区间[a,b上n+1个不同点x0,x1…,xn处的函数值y;=f(x1)(i=0,1,…,n),求一个 至多n次多项式 qn(x)=ao+a1x+…+anx 使其在给定点处与∫(x)同值,即满足插值条件 qn(x)=f(x)=y1(=0,1…,n) φn(x)称为插值多项式,x(=0,1,…n)称为插值节点,简称节点,[a,b]称为插值区 间。从几何上看,n次多项式插值就是过n+1个点(x,f(x,)(i=0,1,…,n),作一条 多项式曲线y=n(x)近似曲线y=f(x)。 n次多项式(1)有n+1个待定系数,由插值条件(2)恰好给出n+1个方程 do +axo taro t.+anxo =yo 10+a1x1+a2x anI=VI ao+a,r, +a2x2+ .+a,rm=y 记此方程组的系数矩阵为A,则 det(a) 是范德蒙特( Vandermonde)行列式。当x0,x1…,xn互不相同时,此行列式值不为零。因 此方程组(3)有唯一解。这表明,只要n+1个节点互不相同,满足插值要求(2)的 插值多项式(1)是唯一的 插值多项式与被插函数之间的差 R, (x)=f(x)-p (x)
-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
称为截断误差,又称为插值余项。当∫(x)充分光滑时, R, (x)=f(x)-L,(x) (n+1))On1(x),s∈(a,b) 其中On1(x)=1(x-x,) 1.12拉格朗日插值多项式 实际上比较方便的作法不是解方程(3)求待定系数,而是先构造一组基函数 (x2-x0)…(x1-x11)(x1-x1+1)…(x1-xn) (i=0,1,…,n) l,(x)是n次多项式,满足 1(x)= L(3)=()=∑ (4) j=0 上式称为n次 Lagrange插值多项式,由方程(3)解的唯一性,n+1个节点的n次 Lagrange插值多项式存在唯一。 1.13用 Matlab作 Lagrange插值 Matlab中没有现成的 Lagrange插值函数,必须编写一个M文件实现 Lagrange插值。 设n个节点数据以数组x0,y0输入(注意 Matlat的数组下标从1开始),m个插值 点以数组x输入,输出数组y为m个插值。编写一个名为 lagrange. n的M文件: function y=lagrange(x0, y0, x)i n=length (x0)im=length(x)i z=x(i) for k=l: n for j=l: n p=p*(z-x0(j))/(x0(k)-x0(j)) end s=p*yo(k)+si end
-176- 称为截断误差，又称为插值余项。当 f (x)充分光滑时， ( ), ( , ) ( 1)! ( ) ( ) ( ) ( ) 1 ( 1) x a b n f R x f x L x n n n n ∈ + = − = + + ω ξ ξ 其中 ∏= + = − n j n j x x x 0 1 ω ( ) ( )。 1.1.2 拉格朗日插值多项式 实际上比较方便的作法不是解方程（3）求待定系数，而是先构造一组基函数 ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) 0 1 1 0 1 1 i i i i i i n i i n i x x x x x x x x x x x x x x x x l x − − − − − − − − = − + − + L L L L , ( 01 ) 0 i , , ,n x x x x n j i j i j j = L − − =∏ ≠ = l (x) i 是n 次多项式，满足 ⎩ ⎨ ⎧ = ≠ = j i j i l x i j 1 0 ( ) 令 ∑ ∑∏ = = ≠ = ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ − − = = n i n i n j i j i j j n i i i x x x x L x y l x y 0 00 ( ) ( ) （4） 上式称为n 次 Lagrange 插值多项式，由方程（3）解的唯一性，n +1个节点的n 次 Lagrange 插值多项式存在唯一。 1.1.3 用 Matlab 作 Lagrange 插值 Matlab 中没有现成的 Lagrange 插值函数，必须编写一个 M 文件实现 Lagrange 插值。 设n 个节点数据以数组 x0, y0输入（注意 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
插值 在导出 Newton公式前,先介绍公式表示中所需要用到的差商、差分的概念及性质。 12.1差商 定义设有函数∫(x),x0,x1 为一系列互不相等的点,称 f(x)-f(x (≠为f(x)关于点x,x,一阶差商(也称均差)记为∫[x,x,] f(x)-f(x,) fIx,, 称一阶差商的差商 fIx x]-f[x, xk] 为f(x)关于点x,x,x的二阶差商,记为fx,x,x]。一般地,称 为f(x)关于点x0,x1…,xk的k阶差商,记为 ∫xnx…J[x,x1…,x--fx2x…,x Ro-x 容易证明,差商具有下述性质 fi,x]=fIxi xi fx,x xk]=f[i,xk, x]=fIx,x,xgI 122 Newton插值公式 线性插值公式可表成 q1(x)=f(x0)+(x-x0)f[x0,x1] 称为一次 Newton插值多项式。一般地,由各阶差商的定义,依次可得 f(x)=f(o) x,Io f[x,x]=f[x0,x1]+(x-x1)[x,x,x1] f[x,x0,x1]=f[x,x1,x2]+(x-x2)f[x,x0,x12x2 ∫[x,x…,xn]=f[x,x,…,xn]+(x-xn)f[x,x0,…,xn] 将以上各式分别乘以1,(x-x0),(x-x0)(x-x1)…,(x-x0)(x-x1)…(x-xm1),然 后相加并消去两边相等的部分,即得 f(x)=f(x0)+(x-x0)f[x0,x1]+ +(x-x0(x-x1)…(x-xn1)[x0,x1,…,xn] +(x-x0(x-x1)…(x-xn)f[x,x0,x1
-177- 1.2 牛顿（Newton）插值 在导出 Newton 公式前，先介绍公式表示中所需要用到的差商、差分的概念及性质。 1.2.1 差商 定义 设有函数 f (x), x0 , x1, x2 ,L为一系列互不相等的点，称 i j i j x x f x f x − ( ) − ( ) (i ≠ j) 为 f (x)关于点 i j x , x 一阶差商（也称均差）记为 [ , ] i j f x x ，即 i j i j i j x x f x f x f x x − − = ( ) ( ) [ , ] 称一阶差商的差商 i k i j j k x x f x x f x x − [ , ]− [ , ] 为 f (x)关于点 i j k x , x , x 的二阶差商，记为 [ , , ] i j k f x x x 。一般地，称 k k k x x f x x x f x x x − − − 0 0 1 1 1 2 [ , ,L, ] [ , ,L, ] 为 f (x)关于点 k x , x , , x 0 1 L 的k 阶差商，记为 k k k k x x f x x x f x x x f x x x − − = − 0 0 1 1 1 2 0 1 [ , , , ] [ , , , ] [ , , , ] L L L 容易证明，差商具有下述性质： [ , ] [ , ] i j j i f x x = f x x [ , , ] [ , , ] [ , , ] i j k i k j j i k f x x x = f x x x = f x x x 1.2.2 Newton 插值公式 线性插值公式可表成 ( ) ( ) ( ) [ , ] 1 0 0 0 1 ϕ x = f x + x − x f x x 称为一次 Newton 插值多项式。一般地，由各阶差商的定义，依次可得 ( ) ( ) ( ) [ , ] 0 0 0 f x = f x + x − x f x x [ , ] [ , ] ( ) [ , , ] 0 0 1 1 0 1 f x x = f x x + x − x f x x x [ , , ] [ , , ] ( ) [ , , , ] 0 1 0 1 2 2 0 1 2 f x x x = f x x x + x − x f x x x x LL [ , , , ] [ , , , ] ( ) [ , , , ] 0 n 1 0 1 n n 0 n f x x L x = f x x L x + x − x f x x L x − 将以上各式分别乘以1,( ),( )( ), , x − x0 x − x0 x − x1 L ( )( ) ( ) − 0 − 1 − n−1 x x x x L x x ，然 后相加并消去两边相等的部分，即得 ( )( ) ( ) [ , , , , ] ( )( ) ( ) [ , , , ] ( ) ( ) ( ) [ , ] 0 1 0 1 0 1 1 0 1 0 0 0 1 n n n n x x x x x x f x x x x x x x x x x f x x x f x f x x x f x x L L L L L + − − − + − − − = + − + − 记
N(x)=f(x0)+(x-x0)[x0,x]+ +(x-x0)(x-x1)…(x-xn-1)f[x0,x1,…,xn Rn(x)=(x-x)(x-x)…(x-xn)[x,x0,x1…,xn] On+1(x)[x,x0,x1…,xn] 显然,N2(x)是至多n次的多项式,且满足插值条件,因而它是f(x)的n次插值 多项式。这种形式的插值多项式称为 Newton插值多项式。Rn(x)称为 Newton插值余 项 Newton插值的优点是:每增加一个节点,插值多项式只增加一项,即 Nn1(x)=Nn(x)+(x-x0)…(x-xn)[x0,x1,…,xn] 因而便于递推运算。而且 Newton插值的计算量小于 Lagrange插值。 由插值多项式的唯一性可知, Newton插值余项与 Lagrange余项也是相等的,即 R (x)=Om+(x)fLx,xo, xu, .,x,] (n+1)!““ )5∈(a,b) 由此可得差商与导数的关系 f(s ni 其中5∈(a,B),a=min{x1},B=max{x;}。 123差分 当节点等距时,即相邻两个节点之差(称为步长)为常数, Newton插值公式的形 式会更简单。此时关于节点间函数的平均变化率(差商)可用函数值之差(差分)来表 定义设有等距节点x=x0+bh(k=0,1,…,n),步长h为常数,f=f(xk)。 称相邻两个节点x,xk+处的函数值的增量f+-f(k=0,1,…,n-1)为函数f(x)在 点x处以h为步长的一阶差分,记为AA,即 Afk=fk1-fK 类似地,定义差分的差分为高阶差分。如二阶差分为 △f=△k4-Mk(k=0,1;…,n-2) 一般地,m阶差分为 △"fk=△"fk+1-△"f(k=2,3,…), 上面定义的各阶差分又称为向前差分。常用的差分还有两种: V=fr -K 称为f(x)在xk处以h为步长的向后差分; 8s=fx 称为f(x)在x处以h为步长的中心差分。一般地,m阶向后差分与m阶中心差分公 式为
-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 阶中心差分公 式为
r=Vm-lfk-Vm-fk d"f=o″-f 差分具有以下性质 (i)各阶差分均可表成函数值的线性组合,例如 A"f=∑(-)|.m v"=∑(-1)|f (i)各种差分之间可以互化。向后差分与中心差分化成向前差分的公式如下 fk=△"f d"fk=△ 124等距节点插值公式 如果插值节点是等距的,则插值公式可用差分表示。设已知节点 xk=x+hh(k=0,1,2,…,m),则有 Nn(x)=f(x0)+[x0,x1](x-x0) +…+∫[Lx0,x1…,xn](x-x0x-x1)…(x-xn-1) f △"f0 (x-x0)+……+ (x-x0)(x-x1)…( h nl h 若令x=x+h,则上式又可变形为 N(xo+h)=J+1△o+…+ r(t-1)…(t-n+1) △"f0 n 上式称为 Newton向前插值公式 13分段线性插值 13.1插值多项式的振荡 用 Lagrange插值多项式Ln(x)近似f(x)(a≤x≤b),虽然随着节点个数的增加, Ln(x)的次数n变大,多数情况下误差|Rn(x)会变小。但是n增大时,Ln(x)的光滑 性变坏,有时会出现很大的振荡。理论上,当n→∞,在,b内并不能保证Ln(x)处 处收敛于f(x)。 Runge给出了一个有名的例子 x∈ 对于较大的|x|,随着n的增大,Ln(x)振荡越来越大,事实上可以证明,仅当|x3.63 时,才有 lim L(x)=f(x),而在此区间外,Ln(x)是发散的。 高次插值多项式的这些缺陷,促使人们转而寻求简单的低次多项式插值。 13.2分段线性插值 简单地说,将每两个相邻的节点用直线连起来,如此形成的一条折线就是分段线性 插值函数,记作l(x),它满足n(x1)=y,且l(x)在每个小区间[x,x+1]上是线性
-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 上是线性
函数(i=01,…,n)。 ln(x)可以表示为 ln(x)=∑y(x) x-x1,x∈[x1,x](=0时舍去) I,(x) ,x∈[x1,x1](=n时舍去) 其它 n(x)有良好的收敛性,即对于x∈[a,b]有, lim I(x)=f(x)。 用Ln(x)计算x点的插值时,只用到x左右的两个节点,计算量与节点个数n无关 但n越大,分段越多,插值误差越小。实际上用函数表作插值计算时,分段线性插值就 足够了,如数学、物理中用的特殊函数表,数理统计中用的概率分布表等 1.3.3用 Matlab实现分段线性插值 用 Matlab实现分段线性插值不需要编制函数程序, Matlab中有现成的一维插值函 数 y=interpl (xo, yo, x, 'method,) method指定插值的方法,默认为线性插值。其值可为 nearest'最近项插值 ' linear''线性插值 spline'逐段3次样条插值 cubic'保凹凸性3次插值。 所有的插值方法要求x0是单调的 当x0为等距时可以用快速插值法,使用快速插值法的格式为'* nearest'、"* linea * spline、" cubic'。 4埃尔米特( Hermite)插值 141 Hermite插值多项式 如果对插值函数,不仅要求它在节点处与函数同值,而且要求它与函数有相同的 阶、二阶甚至更高阶的导数值,这就是 Hermite插值问题。本节主要讨论在节点处插值 函数与函数的值及一阶导数值均相等的 Hermite插值。 设已知函数y=f(x)在n+1个互异节点x0,x1…xn上的函数值y1=f(x1) (i=0,1,…,n)和导数值y’;=∫(x)(i=0.1,…,n),要求一个至多2n+1次的多项式 H(x),使得 f(x,)=y1H(x)=y(i=0,1, 满足上述条件的多项式H(x)称为 Hermite插值多项式 Hermite插值多项式为
-180- 函数(i = 0,1,L,n) 。 I (x) n 可以表示为 ∑= = n i n i i I x y l x 0 ( ) ( ) ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ∈ = − − ∈ = − − = + + + − − − 其它 时舍去 时舍去 ， ， （ ） ） 0 [ , ] , [ , ] ( 0 ( ) 1 1 1 1 1 1 x x x i n x x x x x x x i x x x x l x i i i i i i i i i i i I (x) n 有良好的收敛性，即对于 x∈[a,b]有， lim I (x) f (x) n n = →∞ 。 用 I (x) n 计算 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 插值。 设已知函数 y = f (x) 在n +1个互异节点 n x , x , , x 0 1 L 上的函数值 ( ) i i y = f x (i = 0,1,L,n) 和导数值 ' '( ) i i y = f x (i = 0,1,L,n) ，要求一个至多2n +1次的多项式 H(x)，使得 i i H(x ) = y i i H'(x ) = y' (i = 0,1,L,n) 满足上述条件的多项式 H(x)称为 Hermite 插值多项式。 Hermite 插值多项式为
H(x)=>hI(x-x)(2a,y-y)+y 其中h2= ∑ 142用 Matlab实现 Hermite插值 Matlab中没有现成的 Hermite插值函数,必须编写一个M文件实现插值。 设n个节点的数据以数组x0(已知点的横坐标),y0(函数值),yl(导数值) 输入(注意 Matlat的数组下标从1开始),m个插值点以数组x输入,输出数组y为m 个插值。编写一个名为 hermite. m的M文件: function y=hermite(x0, yo, yl, x)i n=length(x0)im=length(x)i for k=l: m yy=0.0 for i=l: n for j=l:n h=h+((x(k)-x0(j))/(x0(i)-x0(j)))^2; a=1/(x0(i)-x0(j) end yy=yy+h*((x0(1)-x(k))*(2*a*y0()-y1(i))+y0(i)) y(k)=yy end 15样条插值 许多工程技术中提出的计算问题对插值函数的光滑性有较高要求,如飞机的机翼外 形,内燃机的进、排气门的凸轮曲线,都要求曲线具有较高的光滑程度,不仅要连续, 而且要有连续的曲率,这就导致了样条插值的产生。 1.5.1样条函数的概念 所谓样条( Spline)本来是工程设计中使用的一种绘图工具,它是富有弹性的细木 条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线(称为样条曲线),并 使连接点处有连续的曲率。 数学上将具有一定光滑性的分段多项式称为样条函数。具体地说,给定区间[a,b 的一个分划 <x=b 如果函数s(x)满足 (i)在每个小区间[x1,x-1](=0,1,…,n-1)上s(x)是k次多项式 (i)s(x)在[a,b]上具有k-1阶连续导数。 则称s(x)为关于分划Δ的k次样条函数,其图形称为k次样条曲线。x0,x1…xn称为 样条节点,x12x2,…,xn1称为内节点,x0,x称为边界点,这类样条函数的全体记做 181
-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 称为边界点，这类样条函数的全体记做
Sp(A,k),称为k次样条函数空间。 显然,折线是一次样条曲线 若S(x)∈S2(△,k),则s(x)是关于分划△的k次多项式样条函数。k次多项式样 条函数的一般形式为 B S (x-x,) 其中a(=0,1…,k)和B(=12,…,n-1)均为任意常数,而 x,) (x-x),x≥x 在实际中最常用的是k=2和3的情况,即为二次样条函数和三次样条函数 二次样条函数:对于[a,b上的分划△:a=x<x1<…<xn=b,则 S2(x) (x-x)+∈S(△,2) 其中(x-x)2= (x-x,),x2x (j=1,2,…,n-1)。 三次样条函数:对于[a,b上的分划A:a=x0<x1<…<xn=b,则 S(x)=a+a1x+2x2+2x2+∑(x-x)∈S2(△3), (6) ≥X 其中(x-x) ,(j=1,2,…,n-1)。 X<x 利用样条函数进行插值,即取插值函数为样条函数,称为样条插值。例如分段线性插值 是一次样条插值。下面我们介绍二次、三次样条插值 1.52二次样条函数插值 首先,我们注意到S2(x)∈Sp(△,2)中含有n+2个特定常数,故应需要n+2个插 值条件,因此,二次样条插值问题可分为两类: 问题(1): 已知插值节点x1和相应的函数值y(=0,1,…,n)以及端点x(或x)处的导数 值y(或yn),求S2(x)∈Sp(4,2)使得 s(x)=(=012…,n) s2(x)=y(或sn(x)=yn) 问题(2) 已知插值节点x和相应的导数值y(=0,1,2,…,m)以及端点x0(或xn)处的函 数值y(或yn),求S2(x)∈S2(△,2)使得 s2(x)=y(i=0,1,2,…,n) (8) s2(x0)=y(或s2(xn)=y)
-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）
事实上,可以证明这两类插值问题都是唯一可解的。 对于问题(1),由条件(7) S2(x0) +-a S2(x1)=a0+a1x1+a2x=y1 s2(x)=a0+ax+a2x2+∑(x-x)2=y,(=2,3,…,m) S2(xo)=a,+a2o=y'o 引入记号X=(an、a12a2,月1…,Bn1)为未知向量,C=(y,y1…,yn,y0)为已 知向量。 0 x1=X1 A=1 0 (r, -x-D) 于是,问题转化为求方程组AX=C的解X=(a0,a,a2,B1…,Bn-)的问题, 即可得到二次样条函数S2(x)的表达式 对于问题(2)的情况类似。 1.53三次样条函数插值 由于S(x)∈S(A,3)中含有n+3个待定系数,故应需要n+3个插值条件,已知 插值节点x1和相应的函数值f(x)=y(=0,1,2,…,n),这里提供了n+1个条件,还 需要2个边界条件。 常用的三次样条函数的边界条件有3种类型: (i)s3(a)=y,s3(b)=yn。由这种边界条件建立的样条插值函数称为f(x)的 完备三次样条插值函数。 特别地,yo=y'n=0时,样条曲线在端点处呈水平状态 如果∫(x)不知道,我们可以要求s3(x)与∫(x)在端点处近似相等。这时以 x,x1,x2,x3为节点作一个三次 Newton插值多项式Nn(x),以xn,xn=1,xn2,xn=3作一 个三次 Newton插值多项式N6(x),要求 s'(a)=N,(a), s (b)=N,(b) 由这种边界条件建立的三次样条称为f(x)的 Lagrange三次样条插值函数。 (i)s"3(a)=y"o,S"3(b)=y"3。特别地yn=y"n=0时,称为自然边界条件
-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 时，称为自然边界条件
(ⅲ)s3(a+0)=s3(b-0),s"3(a+0)=s"3(b-0),(这里要求S3(a+0) S3(b-0)此条件称为周期条件。 154三次样条插值在 Matlab中的实现 在 Matlab中数据点称之为断点。如果三次样条插值没有边界条件,最常用的方法, 就是采用非扭结(not-a-knot)条件。这个条件强迫第1个和第2个三次多项式的三阶 导数相等。对最后一个和倒数第2个三次多项式也做同样地处理 Matlab中三次样条插值也有现成的函数 y=spline(xo, yO, x) p=cape (x0, yo, conds), y=ppva(pp, x) 其中xOy0是已知数据点,x是插值点,y是插值点的函数值。 对于三次样条插值,我们提倡使用函数 cape, cApe的返回值是pp形式,要求插 值点的函数值,必须调用函数 ppva p= cape(xOyO):使用默认的边界条件,即 Lagrange边界条件。 pp= cape(xOy0, conds)中的 conds指定插值的边界条件,其值可为 complete边界为一阶导数,即默认的边界条件 not-a-knot非扭结条件 periodic周期条件 边界为二阶导数,二阶导数的值[0,0 variationa"设置边界的二阶导数值为[0,0 对于一些特殊的边界条件,可以通过 conds的一个1×2矩阵来表示, conds元素的 取值为1,2。此时,使用命令 pp=cape(xo, yo ext, conds) 其中 yo ext=[left,yo, right,这里left表示左边界的取值,nght表示右边界的取值。 conds(i)=j的含义是给定端点i的j阶导数,即 conds的第一个元素表示左边界的条 件,第二个元素表示右边界的条件, conds=[2,1表示左边界是二阶导数,右边界是一阶 导数,对应的值由let和 right给出。 详细情况请使用帮助 help cape 例1机床加工 待加工零件的外形根据工艺要求由一组数据(x,y)给出(在平面情况下),用程控 铣床加工时每一刀只能沿x方向和y方向走非常小的一步,这就需要从已知数据得到加 工所要求的步长很小的(x,y)坐标。 表1中给出的x,y数据位于机翼断面的下轮廓线上,假设需要得到x坐标每改变 0.1时的y坐标。试完成加工所需数据,画出曲线,并求出x=0处的曲线斜率和 13≤x≤15范围内y的最小值。 y01.21.7202.1201.8 1.01.6 要求用 Lagrange、分段线性和三次样条三种插值方法计算。 解编写以下程序: clc, clear x0=[035791112131415]; 184
-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];