c语言埃尔米特插值思路,【数学建模算法】(26)插值和拟合：埃尔米特(Hermite)插值和样条插值...

1.埃尔米特(Hermite)插值

1.1.Hermite插值多项式

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

设已知函数

在

个互异节点

上的函数值

和导数值

，要求一个至多

次的多项式

，使得：

满足上述条件的多项式

称为Hermite多项式。

Hermite 插值多项式为：

其中：

1.2.用Matlab实现Hermite插值

Matlab 中没有现成的 Hermite 插值函数，必须编写一个 M 文件实现插值。

设

个节点的数据以数组

(已知点的横坐标)，

(函数值)，

(导数值)输入(注意 Matlat 的数组下标从 1 开始)，

个插值点以数组

输入，输出数组

为

个插值。编写一个名为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)本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线(称为样条曲线)，并使连接点处有连续的曲率。

数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间

的一个分划：

如果函数

满足：

(1)在每个小区间

上

是

次多项式。

(2)

在

上具有

阶连续导数。

则称

为关于分划

的

次样条函数，其图形称为

次样条函数。

称为样条节点，

称为内节点，

称为边界点，这类样条函数的全体记做

，则

是关于分划

的

次多项式样条函数。

次多项式样条函数的一般形式为：

其中

和

均为任意常数，而：

在实际中最常用的是

或3的情况，即为二次样条函数和三次样条函数。

二次样条函数：

对于

上的分划

，则：

其中：

三次样条函数：

对于

上的分划

，则：

其中：

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

2.2.二次样条函数插值

首先，我们注意到

中含有

个特定常数，故应需要

+ n 个插值条件，因此，二次样条插值问题可分为两类：

问题(1)：

已知插值节点

和相应的函数值

以及端点

(或

)处的导数值

(或

)，求

使得：

问题(2)：

已知插值节点

和相应的导数值

以及端点

(或

)处的函数值

(或

)，求

使得：

事实上，可以证明这两类插值问题都是唯一可解的。

对于问题(1)有：

引入记号

为未知向量，

为已知向量。

于是，问题转化为求方程组

的解

的问题，即可得到二次样条函数

的表达式。

对于问题(2)的情况类似。

2.3.三次样条函数插值

由于

中含有

个待定系数。故应需要

个待定系数，已知插值节点

和相应的函数值

，这里提供了

个条件，还需要2个边界条件。

常用的三次样条函数的边界条件有 3 种类型：

(1)

。由这种边界条件建立的样条插值函数称为

的完备三次样条插值函数。

特别地，

时，样条曲线在端点处呈水平状态。

如果

不知道，我们可以要求

与

在端点处近似相等。这时以

为节点作一个三次Newton插值多项式

，以

作一个三次Newton插值多项式

，要求：

(2)

。特别地

，称为自然边界条件。

(3)

，(这里要求)

此条件称为周期条件。

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 的一个

矩阵来表示，conds 元素的取值为1，2。此时，使用命令：

pp=csape(x0,y0_ext,conds)

其中 y0_ext=[left, y0, right]，这里 left 表示左边界的取值，right 表示右边界的取值。

conds(i)=j 的含义是给定端点

的

阶导数，即 conds 的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶导数，对应的值由 left 和 right 给出。

例：机床加工

待加工零件的外形根据工艺要求由一组数据

给出(在平面情况下)，用程控铣床加工时每一刀只能沿

方向和

方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的

坐标。

下表给出的数据

数据位于机翼断面的下轮廓线上，假设需要得到

坐标每改变0.1时的

坐标。试完成加工所需数据，画出曲线，并求出

处的曲线斜率和

范围内

的最小值。

数据表：

0

3

5

7

9

11

12

13

14

15

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)插值和拟合：拉格朗日插值)