【数学建模算法】(26)插值和拟合:埃尔米特(Hermite)插值和样条插值

1.埃尔米特(Hermite)插值

1.1.Hermite插值多项式

如果对插值函数,不仅要求它在节点处与函数同值,而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值,这就是 Hermite 插值问题。本节主要讨论在节点处插值函数与函数的值及一阶导数值均相等的 Hermite 插值。

设已知函数y=f(x)n+1个互异节点x_{0}, x_{1}, \cdots, x_{n}上的函数值y_{i}=f\left(x_{i}\right)(i=0,1, \cdots, n)和导数值y_{i}^{\prime}=f^{\prime}\left(x_{i}\right) \quad(i=0,1, \cdots, n),要求一个至多2 n+1次的多项式H(x),使得:
H\left(x_{i}\right)=y_{i} \quad H^{\prime}\left(x_{i}\right)=y_{i}^{\prime} \quad(i=0,1, \cdots, n)
满足上述条件的多项式H(x)称为Hermite多项式

Hermite 插值多项式为:
H(x)=\sum_{i=0}^{n} h_{i}\left[\left(x_{i}-x\right)\left(2 a_{i} y_{i}-y_{i}^{\prime}\right)+y_{i}\right]
其中:
h_{i}=\prod_{j=0 \atop j \neq i}^{n}\left(\frac{x-x_{j}}{x_{i}-x_{j}}\right)^{2}, \quad a_{i}=\sum_{j=0 \atop j \neq i}^{n} \frac{1}{x_{i}-x_{j}}

1.2.用Matlab实现Hermite插值

Matlab 中没有现成的 Hermite 插值函数,必须编写一个 M 文件实现插值。
n个节点的数据以数组x0(已知点的横坐标),y0(函数值),y1(导数值)输入(注意 Matlat 的数组下标从 1 开始),m个插值点以数组x输入,输出数组ym个插值。编写一个名为herite.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

2.样条插值

许多工程技术中提出的计算问题对插值函数的光滑性有较高要求,如飞机的机翼外形,内燃机的进、排气门的凸轮曲线,都要求曲线具有较高的光滑程度,不仅要连续,而且要有连续的曲率,这就导致了样条插值的产生。

2.1.样条函数的概念

所谓样条(Spline)本来是工程设计中使用的一种绘图工具,它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线(称为样条曲线),并使连接点处有连续的曲率

数学上将具有一定光滑性的分段多项式称为样条函数。具体地说,给定区间[a, b]的一个分划:
\Delta : \quad a=x_{0}<x_{1}<\cdots<x_{n-1}<x_{n}=b
如果函数s(x)满足:
(1)在每个小区间\left[x_{i}, x_{i-1}\right](i=0,1, \cdots, n-1)s(x)k次多项式。
(2)s(x)[a, b]上具有k-1
阶连续导数。
则称s(x)为关于分划\Deltak次样条函数,其图形称为k次样条函数。x_{0}, x_{1}, \cdots, x_{n}称为样条节点,x_{1}, x_{2}, \cdots, x_{n-1}称为内节点,x_{0}, x_{n}称为边界点,这类样条函数的全体记做S_{P}(\Delta, k),则s(x)是关于分划\Deltak次多项式样条函数。k次多项式样条函数的一般形式为:
s_{k}(x)=\sum_{i=0}^{k} \frac{\alpha_{i} x^{i}}{i !}+\sum_{j=1}^{n-1} \frac{\beta_{j}}{k !}\left(x-x_{j}\right)_{+}^{k}
其中\alpha_{i}(i=0,1, \cdots, k)\beta_{j}(j=1,2, \cdots, n-1)均为任意常数,而:
\left(x-x_{j}\right)_{+}^{k}=\left\{\begin{array}{l}{\left(x-x_{j}\right)^{k}, x \geq x_{j}} \\ {0, \quad x<x_{j}}\end{array}, \quad(j=1,2, \cdots, n-1)\right.
在实际中最常用的是k=2或3的情况,即为二次样条函数和三次样条函数。

二次样条函数:
对于[a, b]上的分划\Delta : a=x_{0}<x_{1}<\cdots<x_{n}=b,则:
s_{2}(x)=\alpha_{0}+\alpha_{1} x+\frac{\alpha_{2}}{2 !} x^{2}+\sum_{j=1}^{n-1} \frac{\beta_{j}}{2 !}\left(x-x_{j}\right)_{+}^{2} \in S_{P}(\Delta, 2)
其中:
\left(x-x_{j}\right)_{+}^{2}=\left\{\begin{array}{l}{\left(x-x_{j}\right)^{2}, x \geq x_{j}} \\ {0, \quad x<x_{j}}\end{array}, \quad(j=1,2, \cdots, n-1)\right.

三次样条函数:
对于[a, b]上的分划\Delta : a=x_{0}<x_{1}<\cdots<x_{n}=b,则:
s_{2}(x)=\alpha_{0}+\alpha_{1} x+\frac{\alpha_{2}}{2 !} x^{2}+\sum_{j=1}^{n-1} \frac{\beta_{j}}{2 !}\left(x-x_{j}\right)_{+}^{2} \in S_{P}(\Delta, 2)
其中:
\left(x-x_{j}\right)_{+}^{3}=\left\{\begin{array}{l}{\left(x-x_{j}\right)^{3}, x \geq x_{j}} \\ {0, \quad x<x_{j}}\end{array}, \quad(j=1,2, \cdots, n-1)\right.

利用样条函数进行插值,即取插值函数为样条函数,称为样条插值。例如分段线性插值是一次样条插值。下面我们介绍二次、三次样条插值。

2.2.二次样条函数插值

首先,我们注意到s_{2}(x) \in S_{P}(\Delta, 2)中含有n+2个特定常数,故应需要n+2 + n 个插值条件,因此,二次样条插值问题可分为两类:

问题(1):
已知插值节点x_{i}和相应的函数值y_{i}(i=0,1, \cdots, n)以及端点x_{0}(或x_{n})处的导数值y_{0}^{\prime}(或y_{n}^{\prime}),求s_{2}(x) \in S_{p}(\Delta, 2)使得:
\left\{\begin{array}{l}{s_{2}\left(x_{i}\right)=y_{i}(i=0,1,2, \cdots, n)} \\ {s_{2}^{\prime}\left(x_{0}\right)=y_{0}\left(或s_{n}^{\prime}\left(x_{n}\right)=y_{n}\right)}\end{array}\right.

问题(2):
已知插值节点x_{i}和相应的导数值y_{i}^{\prime}(i=0,1,2, \cdots, n)以及端点x_{0}(或x_{n})处的函数值y_{0}(或y_{n}),求s_{2}(x) \in S_{p}(\Delta, 2)使得:
\left\{\begin{array}{l}{s_{2}^{\prime}\left(x_{i}\right)=y_{i}^{\prime}(i=0,1,2, \cdots, n)} \\ {s_{2}\left(x_{0}\right)=y_{0}\left(或s_{n}\left(x_{n}\right)=y_{n}\right)}\end{array}\right.

事实上,可以证明这两类插值问题都是唯一可解的。
对于问题(1)有:
\left\{\begin{array}{l}{s_{2}\left(x_{0}\right)=\alpha_{0}+\alpha_{1} x_{0}+\frac{1}{2} \alpha_{2} x_{0}^{2}=y_{0}} \\ {s_{2}\left(x_{1}\right)=\alpha_{0}+\alpha_{1} x_{1}+\frac{1}{2} \alpha_{2} x_{1}^{2}=y_{1}} \\ {s_{2}\left(x_{j}\right)=\alpha_{0}+\alpha_{1} x_{j}+\frac{1}{2} \alpha_{2} x_{j}^{2}+\frac{1}{2} \sum_{i=1}^{k-1} \beta_{i}\left(x_{j}-x_{i}\right)^{2}=y_{j} \quad(j=2,3, \cdots, n)} \\ {s_{2}^{\prime}\left(x_{0}\right)=\alpha_{1}+\alpha_{2} x_{0}=y_{0}^{\prime}}\end{array}\right.

引入记号X=\left(\alpha_{0}, \alpha_{1}, \alpha_{2}, \beta_{1}, \cdots, \beta_{n-1}\right)^{T}为未知向量,C=\left(y_{0}, y_{1}, \cdots, y_{n}, y_{0}^{\prime}\right)为已知向量。

A=\left[\begin{array}{cccccc}{1} & {x_{0}} & {\frac{1}{2} x_{0}^{2}} & {0} & {\cdots} & {0} \\ {1} & {x_{1}} & {\frac{1}{2} x_{1}^{2}} & {0} & {\cdots} & {0} \\ {1} & {x_{2}} & {\frac{1}{2} x_{2}^{2}} & {\frac{1}{2}\left(x_{2}-x_{1}\right)^{2}} & {\cdots} & {0} \\ {\vdots} & {\vdots} & {\vdots} & {0} & {\vdots} & {\vdots} \\ {1} & {x_{n}} & {\frac{1}{2} x_{n}^{2}} & {\frac{1}{2}\left(x_{n}-x_{1}\right)^{2}} & {\cdots} & {\frac{1}{2}\left(x_{n}-x_{n-1}\right)^{2}} \\ {0} & {1} & {x_{0}} & {0} & {\cdots} & {0}\end{array}\right]

于是,问题转化为求方程组AX=C的解X=\left(\alpha_{0}, \alpha_{1}, \alpha_{2}, \beta_{1}, \cdots, \beta_{n-1}\right)^{T}的问题,即可得到二次样条函数s_{2}(x)的表达式。
对于问题(2)的情况类似。

2.3.三次样条函数插值

由于s_{3}(x) \in S_{p}(\Delta, 3)中含有n+3个待定系数。故应需要n+3个待定系数,已知插值节点x_{i}和相应的函数值f\left(x_{i}\right)=y_{i}(i=0,1,2, \cdots, n),这里提供了n+1个条件,还需要2个边界条件。

常用的三次样条函数的边界条件有 3 种类型:
(1)s_{3}^{\prime}(a)=y_{0}^{\prime}, s_{3}^{\prime}(b)=y_{n}^{\prime}。由这种边界条件建立的样条插值函数称为f(x)的完备三次样条插值函数。
特别地,y_{0}^{\prime}=y_{n}^{\prime}=0时,样条曲线在端点处呈水平状态。
如果f^{\prime}(x)不知道,我们可以要求s_{3}^{\prime}(x)f^{\prime}(x)在端点处近似相等。这时以x_{0}, x_{1}, x_{2}, x_{3}为节点作一个三次Newton插值多项式N_{a}(x),以x_{n}, x_{n-1}, x_{n-2}, x_{n-3}作一个三次Newton插值多项式N_{b}(x),要求:
s^{\prime}(a)=N_{a}^{\prime}(a), s^{\prime}(b)=N_{b}^{\prime}(b)

(2)s^{\prime \prime}(a)=y^{\prime \prime}_{0}, s_{3}^{\prime \prime}(b)=y^{\prime \prime}_{3}。特别地y_{n}^{\prime \prime}=y_{n}^{\prime \prime}=0,称为自然边界条件。

(3)s_{3}^{\prime}(a+0)=s_{3}^{\prime}(b-0), s_{3}^{\prime \prime}(a+0)=s_{3}^{\prime \prime}(b-0),(这里要求)s_{3}(a+0)=s_{3}(b-0)此条件称为周期条件。

2.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 指定插值的边界条件,其值可为:

conds 作用
'complete' 边界为一阶导数,即默认的边界条件
'not-a-knot' 非扭结条件
'periodic' 周期条件
'second' 边界为二阶导数,二阶导数的值[0,0]

对于一些特殊的边界条件,可以通过 conds 的一个1 \times 2矩阵来表示,conds 元素的取值为1,2。此时,使用命令:

pp=csape(x0,y0_ext,conds)

其中 y0_ext=[left, y0, right],这里 left 表示左边界的取值,right 表示右边界的取值。
conds(i)=j 的含义是给定端点 ij 阶导数,即 conds 的第一个元素表示左边界的条件,第二个元素表示右边界的条件,conds=[2,1]表示左边界是二阶导数,右边界是一阶导数,对应的值由 left 和 right 给出。

:机床加工
待加工零件的外形根据工艺要求由一组数据(x, y)给出(在平面情况下),用程控铣床加工时每一刀只能沿x方向和y方向走非常小的一步,这就需要从已知数据得到加工所要求的步长很小的(x, y)坐标。
下表给出的数据x, y数据位于机翼断面的下轮廓线上,假设需要得到x坐标每改变0.1时的y坐标。试完成加工所需数据,画出曲线,并求出x=0处的曲线斜率和13 \leq x \leq 15范围内y的最小值。
数据表:

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

利用Matlab编程,使用Lagrange,分段线性和三次样条三种插值方法计算。

clc,clear
x0=[0 3 5 7 9 11 12 13 14 15];
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)]

(Lagrange函数请参见我的简书文章【数学建模算法】(23)插值和拟合:拉格朗日插值)

推荐阅读更多精彩内容