课程设计---Hermite 插值法的程序设计及应用

课程设计说明书
题目：Hermite 插值法的程序设计及应用学生姓名：
学院：
班级：
指导教师：
2012年 1月 5日
摘要
Hermite 插值是数值分析中的一个重要内容，在相同的节点下得到比拉格朗日插值更高次的插值多项式，而且，相应的曲线在部分节点处也更光滑.在我们所学课程中，只给出了当所有节点处一阶导数均已知时的Hermite 插值.但在实际应用中，并不是所有节点处的一阶导数都是已知的.为此，通过查阅文献、学习总结，给出了更具一般性的Hermite 插值公式.已有的Hermite 插值公式成为本文所得结果的一个特例.
本次课程设计，对Hermite 插值法进行了总结,包括Hermit插值法的理论推导，不同情形下的例，以及在解决实际问题中的应用.同时也给出了Hermite插值公式的Matlab算法.
关键词Hermite 插值；Matlab 实现；数值分析
引言 (1)
第一章 Hermite插值 (2)
§1.1 Hermite插值的概念 (2)
§1.2 Hermite插值简单情形 (3)
§1.2.1简单情形解的存在性 (3)
§1.2.2 简单情形解的存在唯一性 (5)
§1.2.3插值余项 (5)
§1.3 Hermite插值其他情形................................ . (5)
第二章 Hermite插值的Matlab实现 (9)
§2.1 导数完全情形Hermite插值的Matlab实现................... ..9 §2.2导数不完全情形Hermite插值的Matlab实现.. (10)
§2.3 Hermite插值在实际问题中的应用 (13)
参考文献 (15)
附录A (16)
附录B (17)
附录C (19)
在实际工作中, 人们得到的一些数据通常是一些不连续的点, 在土木工程、流体力学、经济学和空气动力学等学科中经常要遇到这样的问题. 此时, 这些数据如果不加以处理, 就难以发现其内在的规律性. 如果用户想得到这些分散点外的其他数值, 就必须运用这些已知的点进行插值.因此，对近似公式的构造产生了插值问题.
在实际问题中，两个变量的关系)(x f y =经常要靠实验和观测来获得，而在通常的情况下只能得到)(x f 在有限个点上的值
.,,1,0),(n i x f y i ==
人们希望找到)(x f 的一个近似函数)(x y φ=，使得
i i y x =)(φ，.,,1,0n i = ○
1 此时，)(x f 称为被插值函数，点n i x x x ,,,0 称为插值结点，)(x φ称为插值函数，○
1为插值条件. 常用的插值法有Lagrange 插值、Newton 插值、最近邻插值、Hermite 插值和三次样条插值插值法等. Lagrange 插值在向量X 区域内的插值较准确, 但向量X 区域之外则不太准确.Newton 插值仅适用于等距节点下的牛顿向前(后) 插值. 最近邻插值是最简便的插值, 在这种算法中, 每一个插值输出像素的值就是在输入图像中与其最临近的采样点的值, 当图像中包含像素之间灰度级变化的细微结构时, 最近邻插值法会在图像中产生人工的痕迹. 最近邻插值的特点是简
单、快速, 缺点是误差较大; 三次样条插值一阶和二阶连续可导, 插值曲线光滑, 插值效果比较好, 应用较广Newton 插值和Lagrange 插值虽然构造比较简单,
但都存在插值曲线在节点处有尖点、不光滑、插值多项式在节点处不可导等缺点.
为了保证插值多项式)(x p n 能更好地逼近)(x f , 对)(x p n 增加一些约束条件, 例如要求)(x p n 在某些结点处与)(x f 的微商相等, 这样就产生了切触插值问题.切触插值即为Hermite 插值.它与被插函数一般有更高的密合度.
本课程设计主要对Hermite 插值法进行总结，对其一般情况，特殊情况进行更进一步的学习，尽量实现其在Matlab 及C++上的程序运行.
第一章 Hermite 插值
实际问题中应用较广为Newton 插值和Lagrange 插值，虽然这辆种插值法构造比较简单, 但都存在插值曲线在节点处有尖点、不光滑、插值多项式在节点处不可导等缺点.为了克这些缺点，我们引入了Hermite 插值.
§1.1 Hermite 插值的概念
定义1.1 许多实际插值问题中，为使插值函数能更好地和原来的函数重合，不但要求二者在节点上函数值相等，而且还要求相切，对应的导数值也相等，甚至要求高阶导数也相等.这类插值称作切触插值，或埃尔米特(Hermite)插值.
该定义给出了Hermite 插值的概念，由此得出Hermite 插值的几何意义，如图1.1.
定义1.2 满足上述要求的插值多项式是埃尔米特插值多项式.记为H (x ). 定义1.3 求一个次数不大于1++r n 的代数多项式 H(x) ,满足:
).
(,,2,1),()(.,,2,1),()(n r r i x f x H n i x f x H i i i i ≤='='== （1-1） 则（1-1）为Hermite 插值条件.
定义1.4 令 ),(22y x ),(33y x ),(44y x
),(11y x
)
,(00y x x
y
图1.1 Hermite 插值多项式的几何意义含义
.)()()()()(00∑∑=='+=r
k k k n k k k x f x x f x x H βα (1-2)
其中，),,1,0)(x (),,1,0)((k n k n k x k ==βα和都是1++r n 次待定多项式并且它们满足如下条件：
⎩
⎨⎧=01)(i k x α k i k i ≠= .,,1,0,n k i = .,,1,0,,,1,0,0)('r i n k x i k ===α
⎩
⎨⎧='01)(i k x β k i k i ≠= .,,1,0,r k i = .,,1,0,,,1,0,0)(n i r k x i k ===β
称（1-2）为Hermite 插值公式.
解决Hermite 插值问题，就是在给定结点处函数值与导数值的基础上根据插值公式构造Hermite 插值多项式，并根据已知条件解出多项式系数.
§1.2 Hermite 插值简单情形
已知函数表： x
0x 1x 2x … m x … n x )(x f
0y 1y 2y … m y … n y )(x f ' 0'y 1'y 2'y … m y ' … n y '
求一个插值多项式,使其满足条件数表.由于数表中包含22+n 个条件,所以能够确定次数不大于12+n 的代数多项式 )(12x H n +.
此情形为导数个数与函数值个数相等的情形,即 Hermite 插值问题的最简单也是最常用情形.
1.2.1简单情形解的存在性
由于Hermite 插值公式（1-2）已给出，接下来只需构造出)(x k α及)(x k β，即认为其存在.在此简介Lagrange-Hermite 插值法构造插值多项式.
Step1 构造)(x k α（n k ,,1,0 =）
由条件)(0)(')(k i x x i k i k ≠==αα知),,,1,0(k i r i x i ≠= 是)(x k α的二重零点.已知Lagrange 插值基函数)(x l k 是n 次多项式，且具有性质
⎩
⎨⎧=≠==i k i k x l ki i k ,1,0)(δ， 则2n 次多项式[]2)(x k k 也具有性质[]ki i k x l δ=2)(，而[]2
)(x l k 的一阶导数在)(k i x i ≠处的值[]()0)()(2)(2='='i k i k i k x l x l x l 所以当k i ≠时，i x 也都是[]2
)(x k k
的两重零点.注意到)(x h k 是12+n 次多项式，而[]2
)(x l k 是n 2次多项式，因此可设),,2,1,0)(()()(2n k x l b ax x k k =+=α其中b a ,为待定常数.显然k i ≠时满足0)(')(==i k i k x x αα，现只要求出b a ,满足k i =时，满足0)(',1)(==k k k k x x αα即可.由此得到确定b a ,的两个方程：
)(2)())(()(2)(1
)()()()(22=+'=++'='=+=+=a x l x al b ax x l x l x b ax x l b ax x k k k k k k k k k k k k k k k k k αα
解出 k k k
k k x x l b x l a ⋅'+='-=)(21)(2 于是[])())((21)(2x l x x x l x k k k k
k -'-=α. Step2 构造)(x k β ),,1,0(n k =
由条件)(0)(')(k i x x i k i k ≠==ββ知),,,1,0(k i r i x i ≠= 是)(x k β的二重零点.因此可设)(x k β也含因子)(2x l k ，又0)(=k k x β，所以)(x k β还含有因式)(k x x -,因此设)()()(2x l x x A x k k k -=β，其中A 为待定常数.
显然)(x k β是12+n 次多项式，且当k i ≠时满足0)(')(==i k i k x x αα，由,1)(='k k
x β可确定A 如下： 1)()(2)()()(2=='⋅⋅-+='A x l x l x x A x Al x k k
k k k k k k k β
所以 )()()(2x l x x x k k k -=β.
到此为止，Hermite 插值问题的解)(12x H n +为
[],)()()())((21)(202
0k k n
k k k k
n k k k k f x l x x f x l x x x l x H '-+-'-=∑∑== 特别地，当=n 1时，满足113003
113003)(,)(,)(,)(y x H y x H y x H y x H '=''='==的三阶Hermite 插值多项式为
+⎪⎪⎭
⎫ ⎝⎛--⎥⎦⎤⎢⎣⎡'-+⎪⎪⎭⎫ ⎝⎛--+=21010000103)(21)(x x x x y x x y x x x x x H 201
0111101)(21⎪⎪⎭⎫ ⎝⎛--⎥⎦⎤⎢⎣⎡'-+⎪⎪⎭⎫ ⎝⎛--+x x x x y x x y x x x x .
§1.2.2 简单情形解的存在唯一性
为了简便理解，下面用流程图来说明解的存在唯一性.详见附录A.
§1.2.3 插值余项
定理 1.1 设)(x f 在包含1+n 个插值结点的最小区间[b a ,]上22+n 次连续可微，则存在与x 有关的ξ，b a <<ξ，使得
),()!
22()()()(222x w n f x H x f n +=-+ξ 其中∏=-=n
0j )()(j x x x w .由此可得到三阶Hermite 插值多项式的误差为：
,)()(!
4)()()()(212043x x x x f x H x f x R --=-=ξ ξ在0x 与1x 之间.
§1.3 Hermite 插值其他情形
已知函数表：
x 0x
1x … m x … n x y
0y 1y … m y … n y
y ' 0y ' 1y ' … m y '
求一个插值多项式,使其满足条件数表.该问题中，导数个数与函数值个数不相等.我们称之为Hermite 插值中其他情形.在此简介Newton-Hermite 插值法构造插值多项式.
先分析插值条件的个数：2++m n 个，那么，所构造的多项式的次数一般不能超1++m n .于是，按牛顿差值的思想，可设
);())(()(),
()()()(1011n n n m n x x x x x x x x x P x N x H ---=+=++ ωω
其中，)(x N n 为n 次牛顿差值多项式;)(x P m 为待定的次数不超过m 次的多项式. 显然：
n i x f x N x H i i n i ,,2,1,0),()()( ===
为确定)(x P m ，对)(x H 求导：
)()()()()()(11x x P x x P x N x H n m n m n
++'+'+'='ωω 根据插值条件)()(i i x f x H '='，有
)()()()()()()()()(111i n i m i n
i n i m i n i m i n i x x P x N x x P x x P x N x H +++'+'='+'+'='ωωω 得到
m i x x N x f x P i n
i n i i m ,,2,1,0,)()()()(1 =''-'=+ω 于是，把求)(x P m 的问题转化为又一个插值问题
已知)(x P m 的函数表 x
1x 2x … m x )(x P m )(1x P m )(2x P m … )(m m x P
确定一个次数不超过m 的插值多项式)(x L m ，使其满足)()(i m i m x P x L =. 根据牛顿差值公式.
)())(](,,[)](,[)()(10000100----++-+=m m m m m m x x x x x x x x P x x x x P x P x P
将上式带回，即得到满足条件
;,,2,1,0),()(;
,,2,1,0),()(m k x f x H n k x f x H k k k k ='='==
的Newton-Hermite 插值多项式.
例1.1 已知函数表： x 0x
1x y 0y
1y y ' 0'y
求一个插值多项式H (x )，使其满足条件：
),()(),()(),()(001100x f x H x f x H x f x H '='==
该问题中，导数个数与函数值个数不相等.我们称之为Hermite 插值中其他情形.在此简介Newton-Hermite 插值法构造插值多项式.
先由函数表
x
x 0 x 1 y
y 0 y 1
作线性插值，即为 []
()01001,)()(x x x x f x f x P -+= 再注意到H (x )与P 1 (x )在节点x 0, x 1上函数值相同，即：
1
1110
010)()()()(y x P x H y x P x H ====
于是，它们的差可以设为 ))(()()(101x x x x K x P x H --=-
其中K 为待定常数，上式又可记为：
))(()()(101x x x x K x P x H --+= (1-3)
为确定K ,对上式求导：
)()()(101x x x x K x P x H -+-+'='
令x = x 0，代入上式，并且注意到插值条件0
0)(y x H '='得： []010*******)(,)()()(y x x K x x f x x K x P x H '=-+=-+'='
于是有
[]0
10
10x x y x x f K -'--=
将上式代入（1-3）得
[]))(()()(100
10
101x x x x x x y x x f x P x H ---'--+
=
[][]))(()(,)(100
10
100100x x x x x x y x x f x x x x f x f ---'--+
-+= (1-4)
可以验证（1-4）所确定的H (x )确实满足插值条件（1-1）.同时也可以看到，构造牛顿——埃米尔特插值多项式，完全采用牛顿插值的构造思想.
最后，也可以把（1-4）式整理成拉格朗日形式：
1001
112
01
0001
10
10
10)()(y x x x
x x x y x
x x x y x
x x x x x x x x x x H '-⎪⎪⎭
⎫ ⎝⎛--+⎪⎪⎭⎫ ⎝⎛--+⎪⎪⎭⎫ ⎝⎛----+-=
插值余项为
()()120)3(2!
3)
()(x x x x f X R --=
ξ， ξ在0x 与1x 之间.
第二章 Hermite 插值的Matlab 实现
§2.1 导数完全情形Hermite 插值的Matlab 实现
在实际应用中，应用最广也是最简单的Hermite 插值情形即为导数完全的情况下，Hermite 插值多项式的拟合.我们首先讨论该情形下的Matlab 程序.
在给出程序之前，我们首先给出该公式所应用的Hermite 插值公式. 定理2.1 设在节点b x x x a n ≤<<≤≤ 21上，
,
)(,)(j j j j y x f y x f '='=,
其中n j ≤≤1，则函数)(x f 在结点处n x x x ,,,21 处的Hermite 插值多项式为
∑=+--=n
i i i i i i i y y y a x x h x y 1])2)([()(
其中 ∑
∏≠=≠=-=--=n
i
j j j
i i n
i
j j j
i j i x x a x x x x h 12
11
;)(
.
该定理的证明详见文献.
该情形下对应的Matlab 程序及流程图详见附录B . 为验证该程序的正确性与有效性，下面给出例2.1. 例2.1 设有如下数据表:
x
0 0.5 1 1.5 2 2.5 3 3.5
)sin(x y = 0 0.4794 0.8145 0.9975 0.9093 0.5985 0.1411 -0.3508 )cos(x y =' 1 0.8776 0.5403 0.0707 -0.4161 -0.8011 -0.9900 -0.9365
在Matlab 工作台输入如下命令：
>> x0=[0,0.5,1,1.5,2,2.5,3,3.5];
y0=[0,0.4794,0.8415,0.9975,0.9093,0.5985,0.1411,- 0.3508]; y1=[1,0.8776,0.5403,0.0707,-0.4161,-0.8011,-0.9900,-0.9365]; x=x0;
y=hermite(x0,y0,y1,x); y
plot(x,y) y2=sin(x); hold on
plot(x,y2,'*r') 则输出结点处的插值：
y =0 0.4794 0.8415 0.9975 0.9093 0.5985 0.1411 -0.3508
)sin(x y =的Hermite 插值多项式的拟合图像如图：
§2.2导数不完全情形Hermite 插值的Matlab 实现
在实际应用中，并不是所有节点处的一阶导数都是已知的，为此，我们给出了更具一般性的Hermite 插值公式及其算法实现，已有的Hermite 插值公式成为本文所得结果的一个特例.
在此首先给出求解Hermite 插值问题的一般性公式。

定理2.1 设在节点b x x x a n ≤<<≤≤ 21上，
}),,1{,,,1(,)()
,,1(,)(n i m k y x f n j y x f k k
i j j k ∉='='==,
图2.1 )sin(x y =的Hermite 插值多项式拟合图像
其中n m ≤≤1，则函数)(x f 在结点处n x x x ,,,21 处的Hermite 插值多项式为
.)()()()()(0
∑∑=='+=r
k k k n
k i i x f x x f x x H βα
其中
{}(){}⎪⎪⎪⎩⎪
⎪⎪⎨⎧
∈⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-⎪⎪⎪⎭⎫ ⎝⎛-+--⋅--⋅--∉--⋅--=∏∑∑∏∏∏≠=≠=≠=≠=≠==m i i k m i m i i k i i n i j j j i i i i n i
j j j i j m n i
j j m k i i i j i j i k k k k k k k i i i i x x x x x
x x x x x x x x x i i i i x x x x x x x x x 1211112111,,,,111,,,,,)( α .,,1n i =
.,,1),
()(11
m k x x x x x x x x x x x k l
k l k i m
k
l l i i i n
i
j j j
i j k =-⋅--⋅--=∏
∏
≠=≠=β
该定理的证明详见文献.
由该公式得到的更具一般性的Hermite 插值公式Matlab 程序详见附录C. 为了验证该程序的一般性，首先给出导数完全情况下的例.
例2.2 已知数表：
x 1 1.2 1.4 1.6 1.8 y
1 1.0954 1.183
2 1.2649 1.3416 y '
0.5
0.4564 0.4226 0.3953 0.3727
根据所列的数据点求出其Hermite 插值多项式，并计算当x=2.0 时的y 值.
在Matlab 中键入下面命令
x=1:0.2:1.8; y=[1, 1.0954, 1.1832, 1.2649, 1.3416]; x1=x; y1=[0.5, 0.4564,0.4226,0.3953,0.3727]; [h,yy]=HermiteInt1(x,y,x1,y1,2)
按回车键得到插值多项式：
h =
43215754781469129/1099511627776000-1145972591322841157/4398046511104000*t+163294776469783529843/211106232532992000*t^2+1665485238488168375/60798594969501696*t^8+15318604908211965505/30399297484750848*t^6-469
1451189851556625/30399297484750848*t^7-50392743046091368807/37999121855938560*t^3-63885352929874938617/60798594969501696*t^5+441043429159790924983/303992974847508480*t^4-130567005823358125/60798594969501696*t^9 yy =
1.4112
该例所得插值多项式的对应的函数图如下图所示：
1.2
1.4 1.6 1.82
2.2
0.6
0.8
1.2
1.4
显然，在导数完全的情况下该程序能够准确的求得Hermite 插值多项.为了体现其一般性，给出例2.3.
例2.3 找出次数小于4的)(x p 满足数表
x 0 1 2 y
0 1 1 y
1
的Hermite 插值多项式.
在Matlab 的命令窗口直接输入以下命令: x=[0,1,2];y=[0,1,1];x1=[0,1];y1=[0,1]; h=HermiteInt1(x,y,x1,y1)
按回车键后得到插值多项式为: h =1/4*t^4-3/2*t^3+9/4*t^2.
该例所得插值多项式对应的函数图如下：
图2.2 例2.2的插值拟合图像
1234
0.5
11.522.53
§2.3 Hermite 插值在实际问题中的应用
实际问题 技术员为了描绘出车门的曲线，他得到数据
x
0 1 2 3 4 5 6 7 8 9 10 y
2.51
3.30
4.04 4.70
5.22 5.54 5.78 5.40 5.57 5.70 5.80 y '
0.8 0.2
分析他如何用数值方法实现.
首先用2.2节所得程序求得插值多项式。

在Matlb 中输入以下命令：
x=[0,10,1,2,3,4,5,6,7,8,9];
y=[2.51,5.80,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70]; x1=[0,10]; y1=[0.8,0.2];
h=HermiteInt1(x,y,x1,y1)
按回车键可得插值多项式为：
h=4/5*t-1132854469/317520000*t^2+251/100+4337490841/423360000*t^3-1714571473981/142884000000*t^4+10955247103/1428840000*t^5-57148609549/19051200000*t^6+33010812389/43545600000*t^7-192428418647/1524096000000*t^8+419662093/30481920000*t^9+33910721/914457600000*t^11-432289699/457228800000*t^10-965039/1524096000000*t^12
为验证该Hermite 插值多项式的拟合程度下面将=t 1，2,, 10将代入该多项式可得到拟合的y 值.
y =2.5100,3.3000,4.0400,4.7000,5.2200,5.5400,5.7800,5.4000,5.5700, 5.7000, 5.8000
图2.3 例2.3的插值拟合图
因此可以看出它与函数表完全拟合.
由插值多项式得到的车门曲线如图2.4与图2.5所示.
1
234
2.5
3.5
44.5
55.5
2
468
4.6
4.8
5.2
5.45.65.8
图2.4 在区间[0，5]的车门曲线图
图2.5 在区间[0，10]的车门曲线图
参考文献
[1]李庆扬，王能超，易大义.数值分析.4 版［M］.北京：清华大学出版社，2001:41-45.
[2]应玮婷，王洁.Hermite插值公式的推广及其matlab 实现.[A].浙江临海:台州学院数学与信息工程学院,2010:03-0004-06
[3]黄明游，刘播，徐涛，数值计算方法[M]，北京：科学出版社.2005.
假设存在两个Hermite 多项式1212,++n n H H
)()()(1212i n i i n x H x f x H ++==
.,,1,0n i =
求导
)(')(')('1212i n i i n x H x f x H ++==
.,,1,0n i =
令
)
()();()()(121212x P x F x H x H x F n n n +++∈-=
.
0)(';
0)(==x F x F
又由+代数定理知)(x F 最多有12+n 个根 n x x x ,,,10 为)(x F 的二重根，即)(x F 至少有22+n 个根
矛盾
图1 解的存在唯一性证明流程图
导数完全的情况下Hermite插值多项式拟合Matlab程序: function f = Hermite(x,y,y_1,x0)
syms t;
f = 0.0;
if(length(x) == length(y))
if(length(y) == length(y_1))
n = length(x);
else
disp('y和y的导数的维数不相等！');
return;
end
else
disp('x和y的维数不相等！');
return;
end
for i=1:n
h = 1.0;
a = 0.0;
for j=1:n
if( j ~= i)
h = h*(t-x(j))^2/((x(i)-x(j))^2);
a = a + 1/(x(i)-x(j));
end
end
f = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));
if(i==n)
if(nargin == 4)
f = subs(f,'t',x0);
else
f = vpa(f,6);
end
end
end
对应流程程图为：
开始
输入100,,y y x 基所需求的节点的值
y x ,维数相等
y y ,维数相
根据Hermite 插值公式计算i h 与i a
得出Hermite 插值多项式
输出x 结点处插值多项式的值
结束
yes
yes
no
no 图2 导数完全的情况Hermite 插值Matlab 程序
附录C
更具一般性的Hermite插值公式的Matlab程序：
function[h,yy]=HermiteInt1(x,y,x1,y1,xx)
%求Hermite 插值.x 为插值节点，y 为相应的函数值；在节点x1 的一阶导数为y1；xx 为插值点.
%输出Hermite 插值函数的表达式h，若输入参数中有插值点xx 时，再输出xx 相应的插值函数值yy.
n=length(x);m=length(x1);
syms t
yy=0;
for i=1:n
%下面求y（i）前的系数
I=0;
%下面这个循环是要找出x（i）在数组x1 中的位置
for j=1:m
if x(i)==x(j)
I=j;
break;
end
end
l1=1;l2=0;
for j=1:n
if j~=i
l1=l1*(t-x(j))/(x(i)-x(j));
l2=l2+1/(x(i)-x(j));
end
end
for j=1:m
if j~=I
l1=l1*(t-x1(j))/(x(i)-x1(j));
l2=l2+1/(x(i)-x1(j));
end
end
if I==0
l2=0;
end
yy=yy+l1*(-l2*(t-x(i))+1)*y(i);
end
for i=1:m
%下面求y1（i）前的系数
l3=1;
for j=1:n
if x(j)~=x1(i)
l3=l3*(t-x(j))/(x1(i)-x(j)); end
end
for j=1:m
if x1(j)~=x1(i)
l3=l3*(t-x1(j))/(x1(i)-x1(j)); end
end
yy=yy+l3*(t-x1(i))*y1(i);
end
h=simplify(yy);
if nargin==5
yy=eval(subs(h,'xx','t'));
end。