数值计算(数值分析)实验4-分段三次埃尔米特(hermite)插值【c程序实现+流程图】
课程设计---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 xy图1.1 Hermite 插值多项式的几何意义含义.)()()()()(00∑∑=='+=rk 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 插值简单情形已知函数表: x0x 1x 2x … m x … n x )(x f0y 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 kk k x x l b x l a ⋅'+='-=)(21)(2 于是[])())((21)(2x l x x x l x k k k kk -'-=α. 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 kx β可确定A 如下: 1)()(2)()()(2=='⋅⋅-+='A x l x l x x A x Al x k kk k k k k k k β所以 )()()(2x l x x x k k k -=β.到此为止,Hermite 插值问题的解)(12x H n +为[],)()()())((21)(2020k k nk k k kn k k k k f x l x x f x l x x x l x H '-+-'-=∑∑== 特别地,当=n 1时,满足113003113003)(,)(,)(,)(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 2010111101)(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 +=-+ξ 其中∏=-=n0j )()(j x x x w .由此可得到三阶Hermite 插值多项式的误差为:,)()(!4)()()()(212043x x x x f x H x f x R --=-=ξ ξ在0x 与1x 之间.§1.3 Hermite 插值其他情形已知函数表:x 0x1x … m x … n x y0y 1y … m y … n yy ' 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 ni 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 ni n i i m ,,2,1,0,)()()()(1 =''-'=+ω 于是,把求)(x P m 的问题转化为又一个插值问题已知)(x P m 的函数表 x1x 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 0x1x y 0y1y y ' 0'y求一个插值多项式H (x ),使其满足条件:),()(),()(),()(001100x f x H x f x H x f x H '='==该问题中,导数个数与函数值个数不相等.我们称之为Hermite 插值中其他情形.在此简介Newton-Hermite 插值法构造插值多项式.先由函数表xx 0 x 1 yy 0 y 1作线性插值,即为 []()01001,)()(x x x x f x f x P -+= 再注意到H (x )与P 1 (x )在节点x 0, x 1上函数值相同,即:11110010)()()()(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,代入上式,并且注意到插值条件00)(y x H '='得: []010*******)(,)()()(y x x K x x f x x K x P x H '=-+=-+'='于是有[]01010x x y x x f K -'--=将上式代入(1-3)得[]))(()()(10010101x x x x x x y x x f x P x H ---'--+=[][]))(()(,)(10010100100x 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)式整理成拉格朗日形式:1001112010001101010)()(y x x xx x x y xx x x y xx 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 插值多项式为∑=+--=ni i i i i i i y y y a x x h x y 1])2)([()(其中 ∑∏≠=≠=-=--=nij j ji i nij j ji j i x x a x x x x h 1211;)(.该定理的证明详见文献.该情形下对应的Matlab 程序及流程图详见附录B . 为验证该程序的正确性与有效性,下面给出例2.1. 例2.1 设有如下数据表:x0 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); yplot(x,y) y2=sin(x); hold onplot(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 插值问题的一般性公式。
分段三次埃尔米特插值

分段三次埃尔米特插值分段线性插值函数的导)(I x h 数是间断的,若在节点k x (k =0,1,…,n )上除已知函数值外还给k f 出导数值k k m f ='(k =0,1,…,n ),这样就可构造一个导数连续的分段插值函数)(I x h ,它满足条件:(1).[][]),(,)(I 11b a C b a C x h ∈代表区间一[]b a ,阶导数连续的函数集合. (2)k k h f x =)(I ,'')(I k k h f x -(k =0,1,…n ). (3))(I x h 在每个小区间上是三次[]1,+k k x x 多项式.由两点三次插值多项式可以知道在)(I x h 区间上表达[]1,+k k x x 式为'21111211211)()(21*)()21()()(I k k k k k k k k k k k k k k k k k k h f x x x x x x x x x x x x x x f x x x x x x x x x ---+--+--+--+--=++--++++)(+'1121)(+++---k k kk k f x x x x x x )(.分段三次埃尔米特插值比分段线性插值效果明显改善,但这种插值要求给出节点上的导数值,所要提供的信息太多,其光滑度也不高(只有一阶导数连续),改进这种以克服其缺点就导致三次样条插值的提出.三次样条插值上面讨论的分段插值函数都有一致收敛性,但光滑性比较差,对于像告诉飞机的机翼形线,船体放样等型值线往往要求有二阶光滑度,即有二阶连续导数.早起工程师制图是,把富有弹性的细长木条(所谓样条)用压铁固定在样点上,在其他地方让它自由的弯曲,然后画下长条的曲线,称为样条曲线.样条曲线实际上有分段三次曲线并接而成,在连接点即样点上要求二阶导数连续,从而数学上加以概括就得到数学样条这一概念.三次样条函数定义 若函数[]b a C x S ,)(2∈,并且在每个小区间上是[]1,+j j x x 三次多项式,其中是给定b x x x a n =<<<= 10节点,则称是节点)(x S 0x ,1x ,…,n x 上的三次样条函数.若在节点上j x 给定函数值)(j j x f y =(j =0,1,…,n ),并且成立j j y x S =)( (j =0,1,…,n ),(1.1) 则称为三次)(x S 样条插值函数.由定义知道要求出)(x S ,在每个小区间上要确定[]1,+j j x x 4个待定系数,一共有个小n 区间,所以应该确定4个参数n .根据在上二)(x S []b a ,阶导数的连续性,在节点j x (j =1,2,…,n -1)处应该满足连续性的条件)0()0(+=-j j x S x S ,)0()0(''+=-j j x S x S (1.2))0()0(''''+=-j j x S x S .一共有3n -3个条件,再加上要满)(x S 足插值条件(1.1),共有4n -2个条件,因此还需要2个条件才能确定)(x S .通常可以在区间[]b a ,端点0x a =,n x b =上各加上一个条件(称为边界条件),可根据实际问题要求给定.常见的有以下3种;(1)已知两端的一阶导数值,即'00')(f x S =,;')(n n f x S =. (1.3)(2)两端的二阶导数已知,即''00'')(f x S =,'''')(n n f x S =, (1.4)其特殊情况为0)()(''0''==n x S x S . (1.5)(3)当)(x f 是以n x -0x 为周期的周期函数时,则要求也是)(x S 周期函数.这时边界条件应满足)0()0(0-=+n x S x S ,)0()0('0'-=+n x S x S , )0()0(''0''-=+n x S x S . (1.6)而此时(1.1)中n y y =0.这样确定的样条函数称)(x S 为周期样条函数. 埃尔米特插值不少实际问题的插值问题不但要求在节点上函数值相等,而且还要求对应的导数值也相等,甚至要求高阶导数也相等,满足这种要求插值的多项式就是埃尔米特(Hermi t e )插值多项式.下面只讨论函数值与导数值个数一样的情况.设在节点上b x x x a n ≤<<<≤ 10,)(i i x f y =,)('j j x f m =(j =0,1,…,n ),要求插值多项式)(x H ,满足条件j j y x H =)(,j j m x H =)('(j =0,1,…,n ). (1.1)这里给出了2n +2个条件,可唯一确定一个次数不超过2n +1的多项式)()(12x H x H n =+,其形式为12121012)(++++++=n n n x a x a a x H .如果根据条件(1.1)来确定2n +2个系数0a ,1a ,…,12+n a ,显然非常复杂,因此,我们依旧采用拉格朗日插值多项式的基函数的方法.先求插值基函数)(x j α及)(x j β(j =0,1,…,n ),一共有2n +2个,每一个基函数都是2n +1次多项式,且满足条件⎪⎪⎭⎪⎪⎬⎫====⎩⎨⎧=≠==).,,1,0,()(,0)(;0)(,,1,,0)(''n k j x x x k j k j x jk k j k j k j jk k j δββαδα (1.2)于是满足条件(1.1)的插值多项式可以写成)()(12x H x H n +=用插值基函数表示的形式[]∑=-+=nj j j j j n x m x y x H 012)()()(βα. (1.3)由条件(1.2)可以知道,有k k n y x H =-)(12,kn m x H =+)('12,(k =0,1,…,n ).下面的问题就是求满足条件(1.2)的基函数以)(x j α及)(x j β.所以,我们可以利用拉格朗日插值基函数)(x l j .令)()()(2x l b ax x a j j +=,由条件(1.2)有1)()()(2=+=j j j j j x l b ax x α,[]0)()(2)()()(''=++=j j j j j j j j j x l b ax x al x l x α,整理得⎩⎨⎧=+=+0)(21'j j j x l a b ax . 解出)(2'j j x l a -=,)(21'j j j x l x b +=.由于)())(()()())(()()(110110n j j j j j j n j j j x x x x x x x x x x x x x x x x x l --------=+-+- ,利用两边取对数再求导数,有∑≠=-=njk k kj j jx x x l 0'1)(,所以有)()1)(21()(20x l x x x x x a j njk k kj j j ∑≠=---=. (1.4)同理,可以得到)()()(2x l x x x j j j -=β. (1.5)同时还证明满足条件(1.1)的插值多项式是唯一的.用反证法,假设及都满)(12x H n +)(12x H n +足条件(1.1),所以有)()()(1212x H x H x n n ++-=ϕ在每个节点上均有二重根,即)(x ϕ有2n +2重根.但是是不高)(x ϕ于2n +1次的多项式,所以0)(≡x ϕ.唯一性得到证明.。
课程设计---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 xy图1.1 Hermite 插值多项式的几何意义含义.)()()()()(00∑∑=='+=rk 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 插值简单情形已知函数表: x0x 1x 2x … m x … n x )(x f0y 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 kk k x x l b x l a ⋅'+='-=)(21)(2 于是[])())((21)(2x l x x x l x k k k kk -'-=α. 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 kx β可确定A 如下: 1)()(2)()()(2=='⋅⋅-+='A x l x l x x A x Al x k kk k k k k k k β所以 )()()(2x l x x x k k k -=β.到此为止,Hermite 插值问题的解)(12x H n +为[],)()()())((21)(2020k k nk k k kn k k k k f x l x x f x l x x x l x H '-+-'-=∑∑== 特别地,当=n 1时,满足113003113003)(,)(,)(,)(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 2010111101)(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 +=-+ξ 其中∏=-=n0j )()(j x x x w .由此可得到三阶Hermite 插值多项式的误差为:,)()(!4)()()()(212043x x x x f x H x f x R --=-=ξ ξ在0x 与1x 之间.§1.3 Hermite 插值其他情形已知函数表:x 0x1x … m x … n x y0y 1y … m y … n yy ' 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 ni 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 ni n i i m ,,2,1,0,)()()()(1 =''-'=+ω 于是,把求)(x P m 的问题转化为又一个插值问题已知)(x P m 的函数表 x1x 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 0x1x y 0y1y y ' 0'y求一个插值多项式H (x ),使其满足条件:),()(),()(),()(001100x f x H x f x H x f x H '='==该问题中,导数个数与函数值个数不相等.我们称之为Hermite 插值中其他情形.在此简介Newton-Hermite 插值法构造插值多项式.先由函数表xx 0 x 1 yy 0 y 1作线性插值,即为 []()01001,)()(x x x x f x f x P -+= 再注意到H (x )与P 1 (x )在节点x 0, x 1上函数值相同,即:11110010)()()()(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,代入上式,并且注意到插值条件00)(y x H '='得: []010*******)(,)()()(y x x K x x f x x K x P x H '=-+=-+'='于是有[]01010x x y x x f K -'--=将上式代入(1-3)得[]))(()()(10010101x x x x x x y x x f x P x H ---'--+=[][]))(()(,)(10010100100x 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)式整理成拉格朗日形式:1001112010001101010)()(y x x xx x x y xx x x y xx 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 插值多项式为∑=+--=ni i i i i i i y y y a x x h x y 1])2)([()(其中 ∑∏≠=≠=-=--=nij j ji i nij j ji j i x x a x x x x h 1211;)(.该定理的证明详见文献.该情形下对应的Matlab 程序及流程图详见附录B . 为验证该程序的正确性与有效性,下面给出例2.1. 例2.1 设有如下数据表:x0 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); yplot(x,y) y2=sin(x); hold onplot(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 插值问题的一般性公式。
数值分析Hermite

由(2.3)可设
0 x x x1 a x x0 b ,
2
再由(2.2)可求得
b 1
x1 x0
2
, a
2
x1 x0
பைடு நூலகம்
3
x x1 x x0 0 x 1 2 x1 x0 x0 x1
4!
3. 2n+1 次Hermite 插值多项式
给定n+1个节点和相应的函数值和导数值:
f xi yi , yi f xi mi , i 0,1,, n
则可构造2n+1 次Hermite 插值多项式H x 满足条件:
H 1) x 是不超过2n+1 次多项式; H 2) xi yi , H xi mi , i 0,1,, n
2
2
f , x0 , x1 4!
4
于是有下述定理
定理:设 H3 x 是以 x0 , x1 为插值节点的三次 f x C 3 a, b , f 4 x 在 a, b Hermite 插值多项式, 内存在,其中 a, b 是包含 x0 , x1 的任一区 间,则对任意给定的 x a, b ,总存在依赖 于 x 的点 a, b ,使 4 f 2 2 R3 x f x H 3 x x x0 x x1 .
Hermite插值方法

数值分析实验报告五一、实验目的理解Hermite插值方法,掌握Hermite插值算法设计二、实验内容使用vc++编程,实现该方法,即Hermite插值法三、实验步骤#include <iostream.h>double herm(double x0,double x1,double y0,double y1,double h0,double g0,double g1,double x) {double alp0,alp1,bta0,bta1,t;double s;t=h0*h0;alp0=(x-x1)*(x-x1)*(h0+2*(x-x0))/t/h0;alp1=(x-x0)*(x-x0)*(h0-2*(x-x1))/t/h0;bta0=(x-x0)*(x-x1)*(x-x1)/t;bta1=(x-x1)*(x-x0)*(x-x0)/t;s=y0*alp0+y1*alp1+g0*bta0+g1*bta1;return(s);}void main(){int n=7;double p0;double pn; double aa[8],bb[8],s=0;double xx[8]={0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9};double yy[8]={0.4794,0.6442,0.7833,0.8912,0.9636,0.9975,0.9917,0.9463};double g[8];int i;double a[8],c[8],h[8];cout<<"Please input p0 and pn"<<endl;cin>>p0;cin>>pn;for(i=0;i<=n-1;i++){h[i]=xx[i+1]-xx[i];cout<<"h["<<i<<"]="<<h[i]<<endl;}c[0]=1;g[0]=3*(yy[1]-yy[0])/h[0]-p0*h[0]/2;for( i=1;i<=n-1;i++){a[i]=h[i]/(h[i]+h[i-1]);c[i]=1-a[i];}for(i=1;i<n;i++){cout<<"a["<<i<<"]="<<a[i]<<endl;cout<<"c["<<i<<"]="<<c[i]<<endl;}for( i=1;i<=n-1;i++){g[i]=3*(c[i]*(yy[i+1]-yy[i])/h[i]+a[i]*(yy[i]-yy[i-1])/h[i-1]);}a[n]=1;g[n]=3*(yy[n]-yy[n-1])/h[n-1]+pn*h[n-1]/2;for(i=0;i<=n;i++)cout<<"g["<<i<<"]="<<g[i]<<endl;aa[0]=2;bb[0]=c[0]/aa[0];g[0]=g[0]/aa[0];for(i=1;i<=n-1;i++){aa[i]=2-a[i]*bb[i-1];bb[i]=c[i]/aa[i];g[i]=(g[i]-a[i]*g[i-1])/aa[i];}aa[n]=2-a[n]*bb[n-1];g[n]=(g[n]-a[n]*g[n-1])/aa[n];for(i=n-1;i>=0;i--){g[i]=g[i]-bb[i]*g[i+1];}cout<<endl;for(i=0;i<=n;i++)cout<<"g["<<i<<"]="<<g[i]<<endl;double ss;double c0,c1,d0,d1,g0,g1,h1;double x0;cout<<"Please input interpolation point x0:"<<endl;cin>>x0;if(x0>=0.5 && x0<0.7){c0=xx[0];c1=xx[1];d0=yy[0];d1=yy[1];h1=h[0];g0=g[0];g1=g[1];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=0.7 && x0<0.9){c0=xx[1];c1=xx[2];d0=yy[1];d1=yy[2];h1=h[1];g0=g[1];g1=g[2];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=0.9 && x0<=1.1){c0=xx[2];c1=xx[3];d0=yy[2];d1=yy[3];h1=h[2];g0=g[2];g1=g[3];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=1.1 && x0<=1.3){c0=xx[3];c1=xx[4];d0=yy[3];d1=yy[4];h1=h[3];g0=g[3];g1=g[4];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=1.3 && x0<=1.5){c0=xx[4];c1=xx[5];d0=yy[4];d1=yy[5];h1=h[4];g0=g[4];g1=g[5];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=1.5 && x0<=1.7){c0=xx[5];c1=xx[6];d0=yy[5];d1=yy[6];h1=h[5];g0=g[5];g1=g[6];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=1.7 && x0<=1.9){c0=xx[6];c1=xx[7];d0=yy[6];d1=yy[7];h1=h[6];g0=g[6];g1=g[7];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}elsecout<<"The data error,please input again!"<<endl;}四、运行结果。
埃尔米特(Hermite)插值

实验二埃尔米特(Hermite)插值一、实验目的:1.掌握埃尔米特插值算法原理;2.使用C语言编程实现埃尔米特插值算法。
二、实验准备:阅读《数值分析》2.4节二、实验要求:某人从甲地开车去乙地,每隔一段时间对行车距离和速率进行一次采样,得到在n+1 个采样时刻点t i 的里程s i和速率v i(i=0, 1, ..., n)。
要求编程构造埃尔米特插值多项式H2n+1(t),满足H2n+1(t i)=s i,H'2n+1(t i)=v i,对所有i=0, 1, ..., n成立,并据此计算m个给定时刻的里程和速率。
函数接口定义:void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] );其中N为采样点个数(注意这个N不是公式中的最大下标n,而是等于n+1),采样时刻点t i、里程s i、速率v i分别通过t、s、v传入;m是需要估算的给定时刻的个数,ht传入给定的时刻点,相应计算出的里程和速率应分别存储在hs和hv中。
裁判程序如下:裁判输入数据:20.0 1.00.0 1.00.0 0.050.0 0.2 0.5 0.8 1.030.0 0.5 1.0100.0 170.0 200.030.0 150.0 0.050.0 0.25 0.5 0.75 1.050.0 1.0 2.0 3.0 4.00.0 60.0 160.0 260.0 300.05.0 70.0 100.0 120.0 20.0100.5 1.0 1.5 2.0 2.5 3.0 3.5 3.8 3.95 4.0标准输出数据:0.0000 0.1040 0.5000 0.8960 1.00000.0000 0.9600 1.5000 0.9600 0.0000100.0000 127.9297 170.0000 195.9766 200.000030.0000 165.4688 150.0000 52.9688 0.000030.2222 60.0000 105.9303 160.0000 206.3438 260.0000 307.9764 305.7687 299.9796 300.000062.6024 70.0000 109.0488 100.0000 92.9745 120.0000 41.2374 -44.8421 -16.2783 20.0000#include<stdio.h>#define MAXN 5 /* 最大采样点个数 */#define MAXM 10 /* 最大估算点个数 */void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] ){double l[10],p[10],h1[10],h2[10],x,ll[10],pp[10];int kk;for(kk=0;kk<m;kk++){x=ht[kk];hs[kk]=0;hv[kk]=0;int i;for(i=0;i<N;i++){l[i]=1;ll[i]=1;int j;for(j=0;j<N;j++){if(i!=j){l[i]=l[i]*(x-t[j])/(t[i]-t[j]);}}p[i]=0;pp[i]=0;int k;for(k=0;k<N;k++){if(i!=k){p[i]=p[i]+l[i]/(x-t[k]);pp[i]=pp[i]+ll[i]/(t[i]-t[k]);}}h1[i]=(1-2*pp[i]*(x-t[i]))*l[i]*l[i];h2[i]=(x-t[i])*l[i]*l[i];hs[kk]=hs[kk]+s[i]*h1[i]+v[i]*h2[i];int kkk;for(kkk=0;kkk<N;kkk++){if(x==t[kkk])break;}if(x==t[kkk])hv[kk]=v[kkk];elsehv[kk]=hv[kk]+s[i]*(2*p[i]*l[i]-4*l[i]*p[i]*(x-t[i])*pp[i]-2*pp[i]*l[ i]*l[i])+v[i]*(l[i]*l[i]+2*l[i]*p[i]*(x-t[i]));}}}int main(){int N, m;double t[MAXN], s[MAXN], v[MAXN]; /* 用于构造的数据 */double ht[MAXM], hs[MAXM], hv[MAXM]; /* 用估算的数据 */int i;while ( scanf("%d", &N) != EOF ) {for ( i=0; i<N; i++ )scanf("%lf", &t[i]);for ( i=0; i<N; i++ )scanf("%lf", &s[i]);for ( i=0; i<N; i++ )scanf("%lf", &v[i]);scanf("%d", &m);for ( i=0; i<m; i++ )scanf("%lf", &ht[i]);Hermite_Interpolation( N, t, s, v, m, ht, hs, hv );for ( i=0; i<m; i++ )printf("%.4lf ", hs[i]);printf("\n");for ( i=0; i<m; i++ )printf("%.4lf ", hv[i]);printf("\n\n");}return 0; }。
数值计算方法( 三次样条插值)

u xj hj
分段三次Hermite插值算法
则 v A1 y j 1 A2 y j B1 f j1 B2 f j
算法: 1.输入x j , f j , f j (j 0,1,...,n); 2.计算插值 (1)输入插值点u; (2)对于j 1,2,...,n做 如果u x j 则计算A1 , A2 , B1 , B2 ; v A1 f j 1 A2 f j B1 f j1 B2 f j; 3.输出u , v。
三次样条插值
于是由Taylor展示有 s( x) s( xi ) s( xi )(x xi ) s( xi ) s( xi ) 2 ( x xi ) ( x xi )3 2! 3! M M Mi yi s( xi )(x x j ) i ( x xi ) 2 i 1 ( x xi )3 2! 3!( xi 1 xi )
2M 0 M 1 6 f [ x0 , x0 , x1 ]
三次样条插值
同理(2)式中令i n得 M n 1 2M n 6 f [ xn 1 , xn , xn ] 即有 2M 0 M 1 6 f [ x0 , x0 , x1 ] ) i M i 1 2M i i M i 1 6 f [ xi 1 , xi , xi 1 ] (i 1,2,...,n 1 M 2M 6 f [ x , x , x ] n n 1 n n n 1
三次样条插值
对于待定系数a j , b j , c j .d j j 1,2,...n,即4n个未知系数,
而插值条件为 n 2个,还缺两个,因此须 4 给出两个 条件称为边界条件,有 以下三类: 第一类 已知两端点的一阶导数 s( x0 ) f ( x0 ) m0 s( xn ) f ( xn ) mn
埃尔米特(Hermite)插值

实验二埃尔米特(Hermite)插值一、实验目的:1.掌握埃尔米特插值算法原理;2.使用C语言编程实现埃尔米特插值算法。
二、实验准备:阅读《数值分析》2.4节二、实验要求:某人从甲地开车去乙地,每隔一段时间对行车距离和速率进行一次采样,得到在n+1 个采样时刻点t i 的里程s i和速率v i(i=0, 1, ..., n)。
要求编程构造埃尔米特插值多项式H2n+1(t),满足H2n+1(t i)=s i,H'2n+1(t i)=v i,对所有i=0, 1, ..., n成立,并据此计算m个给定时刻的里程和速率。
函数接口定义:void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] );其中N为采样点个数(注意这个N不是公式中的最大下标n,而是等于n+1),采样时刻点t i、里程s i、速率v i分别通过t、s、v传入;m是需要估算的给定时刻的个数,ht传入给定的时刻点,相应计算出的里程和速率应分别存储在hs和hv中。
裁判程序如下:裁判输入数据:20.0 1.00.0 1.00.0 0.050.0 0.2 0.5 0.8 1.030.0 0.5 1.0100.0 170.0 200.030.0 150.0 0.050.0 0.25 0.5 0.75 1.050.0 1.0 2.0 3.0 4.00.0 60.0 160.0 260.0 300.05.0 70.0 100.0 120.0 20.0100.5 1.0 1.5 2.0 2.5 3.0 3.5 3.8 3.95 4.0标准输出数据:0.0000 0.1040 0.5000 0.8960 1.00000.0000 0.9600 1.5000 0.9600 0.0000100.0000 127.9297 170.0000 195.9766 200.000030.0000 165.4688 150.0000 52.9688 0.000030.2222 60.0000 105.9303 160.0000 206.3438 260.0000 307.9764 305.7687 299.9796 300.000062.6024 70.0000 109.0488 100.0000 92.9745 120.0000 41.2374 -44.8421 -16.2783 20.0000#include<stdio.h>#define MAXN 5 /* 最大采样点个数 */#define MAXM 10 /* 最大估算点个数 */void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] ){double l[10],p[10],h1[10],h2[10],x,ll[10],pp[10];int kk;for(kk=0;kk<m;kk++){x=ht[kk];hs[kk]=0;hv[kk]=0;int i;for(i=0;i<N;i++){l[i]=1;ll[i]=1;int j;for(j=0;j<N;j++){if(i!=j){l[i]=l[i]*(x-t[j])/(t[i]-t[j]);}}p[i]=0;pp[i]=0;int k;for(k=0;k<N;k++){if(i!=k){p[i]=p[i]+l[i]/(x-t[k]);pp[i]=pp[i]+ll[i]/(t[i]-t[k]);}}h1[i]=(1-2*pp[i]*(x-t[i]))*l[i]*l[i];h2[i]=(x-t[i])*l[i]*l[i];hs[kk]=hs[kk]+s[i]*h1[i]+v[i]*h2[i];int kkk;for(kkk=0;kkk<N;kkk++){if(x==t[kkk])break;}if(x==t[kkk])hv[kk]=v[kkk];elsehv[kk]=hv[kk]+s[i]*(2*p[i]*l[i]-4*l[i]*p[i]*(x-t[i])*pp[i]-2*pp[i]*l[ i]*l[i])+v[i]*(l[i]*l[i]+2*l[i]*p[i]*(x-t[i]));}}}int main(){int N, m;double t[MAXN], s[MAXN], v[MAXN]; /* 用于构造的数据 */double ht[MAXM], hs[MAXM], hv[MAXM]; /* 用估算的数据 */int i;while ( scanf("%d", &N) != EOF ) {for ( i=0; i<N; i++ )scanf("%lf", &t[i]);for ( i=0; i<N; i++ )scanf("%lf", &s[i]);for ( i=0; i<N; i++ )scanf("%lf", &v[i]);scanf("%d", &m);for ( i=0; i<m; i++ )scanf("%lf", &ht[i]);Hermite_Interpolation( N, t, s, v, m, ht, hs, hv );for ( i=0; i<m; i++ )printf("%.4lf ", hs[i]);printf("\n");for ( i=0; i<m; i++ )printf("%.4lf ", hv[i]);printf("\n\n");}return 0; }。
数值分析实验六(分段三次Hermite插值)

数值分析实验六(分段三次Hermite插值)《数值分析》实验报告实验编号:实验六课题名称:分段三次Hermite插值一、算法介绍给定的函数为f(x)=1/(25*x*x+1),将给定区间分成10分,得到11个节点:x[0],x[1],...,x[10],构造插值函数的基函数。
当x在(x[0],x[1])区间上时,H[0] = (x-x[0])*[((x-x[1])/(x[0]-x[1]))^2]。
其余的区间为H[0]=0。
h[0]= [1+2*(x-x[0])/(x[1]-x[0])]*[((x-x[1])/(x[0]-x[1]))^2]。
当x在[x[i-1],x[i]] (i=1,2,3, (9)区间上时,H[i]=(x-x[i])*[((x-x[i-1])/(x[i]-x[i-1]))^2],h[i]=[1+2*(x-x[i])/(x[i-1]-x[i])]*[((x-x[i-1])/(x[i]-x[i-1]))^2)。
当x在(x[i],x[i+1]](i=1,2,3,…,10)区间上其余的区间为H[i]=(x-x[i])[((x-x[i+1])/(x[i]-x[i+1]))^2],h[i]=[1+2*(x-x[i])/(x[i+1]-x[i])]*[((x-x[i+1 ])/(x[i]-x[i+1]))^2]。
其余区间上均为H[i]=0,h[i]=0(i=1,2,…,10)。
当x在(x[9],x[10])区间上时,H[10] = (x-x[9])(((x-x[10])/(x[9]-x[10]))^2).其余的区间为H[10]=0.h[10]= (1+2*((x-x[9])/(x[10]-x[9])))(((x-x[10])/(x[9]-x[10]))^2).其余区间h[10]=0。
构造函数H(x) =∑(y[i]*h[i]+y'[i]*H[i],(i=0,1,…,10)。
二、程序代码// testV iew.cpp : implementation of the CT estV iew class//#include "stdafx.h"#include "test.h"#include "testDoc.h"#include "testView.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////////////////////////////// //////////// CTestV iewIMPLEMENT_DYNCREA TE(CTestView, CView)BEGIN_MESSAGE_MAP(CTestView, CView)//{{AFX_MSG_MAP(CTestView)// NOTE - the ClassWizard will add and remove mapping macros here.// DO NOT EDIT what you see in these blocks of generated code!//}}AFX_MSG_MAP// Standard printing commandsON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_DIRECT, CV iew::OnFilePrint) ON_COMMAND(ID_FILE_PRINT_PREVIEW,CView::OnFilePrintPreview)END_MESSAGE_MAP()/////////////////////////////////////////////////////////////////// //////////// CTestV iew construction/destructionCTestView::CTestV iew(){// TODO: add construction code here}CTestView::~CT estView(){}BOOL CTestView::PreCreateWindow(CREA TESTRUCT& cs){// TODO: Modify the Window class or styles here by modifying // the CREA TESTRUCT csreturn CV iew::PreCreateWindow(cs);}/////////////////////////////////////////////////////////////////// //////////// CTestV iew drawingvoid CTestView::OnDraw(CDC* pDC){CTestDoc* pDoc = GetDocument();ASSERT_V ALID(pDoc);// TODO: add draw code for native data hereint i,j,k;double x,y,p_x,p_y,l,xx[100],f[100],F[100],sum,p_sum;CPen MyPen,*OldPen;pDC->SetViewportOrg(400,400); //定义坐标原点for(i=-500;i<500;i++){pDC->SetPixel(0,i,RGB(0,0,0));pDC->SetPixel(i,0,RGB(0,0,0)); //画出坐标}pDC->TextOut(-210,5,"-1");pDC->TextOut(196,5,"1");//原函数MyPen.CreatePen(PS_SOLID,1,RGB(255,0,0));//定义画笔颜色OldPen=pDC->SelectObject(&MyPen);x=-1.0,y=1/(1+25*x*x);p_x=x*200;p_y=-y*200;pDC->MoveTo(p_x,p_y);for (x=-1.0;x<=1.0;x+=0.0001){y=1/(1+25*x*x);p_x=x*200;p_y=-y*200;pDC->LineT o(p_x,p_y);}pDC->SelectObject(OldPen);MyPen.DeleteObject();//分段三次Hermite插值MyPen.CreatePen(PS_SOLID,1,RGB(0,0,0)); OldPen=pDC->SelectObject(&MyPen); x=-1.0,y=1.0/(1+25*x*x);p_x=x*200;p_y=-y*200;pDC->MoveTo(p_x,p_y);x=-1.0;for(i=0;i<=10;i++){f[i]=1/(1+25*x*x);xx[i]=x;F[i]=-50*x/(1+25*x*x)/(1+25*x*x); //导数x+=0.2;}x=-1.0;for(j=0;j<=1000;j++){sum=0;for(i=0;i<=10;i++){if(x==xx[i]){sum=f[i];p_x=x*200,p_y=-sum*200;pDC->LineT o(p_x,p_y);break;}if(xxx[i]){y=(1+2*(x-xx[i])/(xx[i+1]-xx[i]))*(x-xx[i+1])*(x-xx[i+1])/(xx[i]-xx[i+1])/(xx[i]-xx[i+1]);sum+=f[i]*y;y=(1+2*(x-xx[i+1])/(xx[i]-xx[i+1]))*(x-xx[i])*(x-xx[i])/(xx[i+1]-xx[i])/(xx[i+1]-xx[i]);sum+=f[i+1]*y;y=(x-xx[i])*(x-xx[i+1])*(x-xx[i+1])/(xx[i]-xx[i+1])/(xx[i]-xx[i+1]);sum+=F[i]*y;y=(x-xx[i+1])*(x-xx[i])*(x-xx[i])/(xx[i+1]-xx[i])/(xx[i+1]-xx[i]);sum+=F[i+1]*y;p_x=x*200;p_y=-sum*200;pDC->LineT o(p_x,p_y);break;}}x+=0.002;}pDC->SelectObject(OldPen);MyPen.DeleteObject();/////////////////////////////////////////////////////////////////// //////////// CTestV iew printingBOOL CTestView::OnPreparePrinting(CPrintInfo* pInfo){// default preparationreturn DoPreparePrinting(pInfo);}void CTestView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add extra initialization before printing}void CTestView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add cleanup after printing}/////////////////////////////////////////////////////////////////// //////////// CTestV iew diagnostics#ifdef _DEBUGvoid CTestView::AssertV alid() const{CView::AssertV alid();}void CTestView::Dump(CDumpContext& dc) const{CView::Dump(dc);CTestDoc* CT estV iew::GetDocument() // non-debug version is inline{ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CT estD oc)));return (CT estDoc*)m_pDocument;}#endif //_DEBUG/////////////////////////////////////////////////////////////////// //////////// CTestV iew message handlers三、运算结果截屏红色的曲线为原函数图像,黑色曲线为分段三次Hermite插值曲线四、算法分析上述图像中黑色的曲线为分段分段三次Hermite插值多项式所对应的图像,由图像可看出黑色的分段三次Hermited插值函数图像和拉格朗日、分段线性插值相比与红色被逼近函数的重合度最好,说明分段三次Hermite插值在函数的各节点两边插值函数的导数是相等的,保证了在各节点的平滑性,且不会出现Runge现象。
三次Hermite插值

04 实例分析
CHAPTER
实例一:已知数据点的插值
总结词
利用已知数据点进行插值,可三次Hermite插值方法,利用已知的数据点来估计未知点的值。这 种方法能够更好地处理数据点的变化,并提高插值的精度。
CHAPTER
插值多项式的构造
定义
Hermite插值法是一种通过已知的离散数据点来构造一个多 项式,使其能够准确地经过这些数据点,并尽可能地平滑地 连接这些点的方法。
构造方法
Hermite插值多项式由两个部分组成,一个是线性函数,另 一个是二次函数。线性函数部分用于确保插值多项式能够准 确地经过数据点,而二次函数部分则用于保证插值多项式的 平滑性。
实例二:未知数据点的插值
总结词
在未知数据点的情况下,可以通过三次 Hermite插值方法,预测并估计未知点的值。
详细描述
在数据点未知的情况下,可以利用三次 Hermite插值方法,根据已知的数据点来预 测和估计未知点的值。这种方法能够为后续 的数据分析和处理提供重要的参考依据。
实例三:复杂函数的插值
三次Hermite插值能够提供高精度的插值结果,特别是在处理
复杂函数时。
稳定性好
02
该方法在处理大数据集时表现出良好的稳定性,不易受到噪声
和异常值的影响。
易于实现
03
三次Hermite插值的算法相对简单,易于在计算机上实现和优
化。
三次Hermite插值的局限性
对初始数据敏感
三次Hermite插值的结果对初始数据的选择 较为敏感,不同的初始数据可能导致不同的 插值结果。
数值计算课后答案4

习题 四 解 答之勘阻及广创作1、设010,1x x ==,写出()x f x e -=的一次插值多项式1()L x ,并估计插值误差。
设插值函数为1()L x ax b =+解之得111a eb -⎧=-⎨=⎩则11()(1)1L x e x -=-+因为(),()x x y x e y x e --'''=-= 所以,插值余项为 所以010101()max max (1)2111248x r x e x x e ξξ-≤≤≤≤-≤-=⨯⨯=。
2解:设三次插值多项式为230123()f x a a x a x a x =+++,由插值条件,建立方程组为 即解之得则所求的三次多项式为23()0.41 6.29 3.489.98f x x x x =--+。
所以3、设(0,1,2,,)i x i n =是 n+1个互异节点,证明: (1)0()(0,1,2,,)nk k i i i x l x x k n ===∑;(2)0()()0(0,1,2,,)n k i i i x x l x k n =-==∑。
证明: (1)由拉格朗日插值定理,以x 0,x 1,x 2,…x n 为插值节点,对y=f(x)=x k作n 次插值,插值多项式为0()()nn i i i p x l x y ==∑,而y i =x i k,所以0()()()nnk n i i i i i i p x l x y l x x ====∑∑同时,插值余项 所以0()nk k i i i l x x x ==∑结论得证。
(2)取函数()(),0,1,2,,k f x x t k n =-=对此函数取节点(0,1,2,,)i x i n =,则对应的插值多项式为0()()()nk n i i i p x x t l x ==-∑,由余项公式,得(1)(1)011()()()()()()()()0(1)!(1)!nn kkn ki i i r x x t x t l x f x x t x n n ξξππ++==---==-=++∑所以令t=x ,4、给定数据(()f x =(1)试用线性插值计算f(2.3)的近似值,并估计误差;(2)试用二次Newton 插值多项式计算f(2.15)的近似值,并估计误差。
数值分析实验报告Hermite插值法、Runge现象,比较Language插值、分段线性插值、分段三次Hermie插值

山东师范大学数学科学学院实验报告x 0.1 0.5 1 1.5 2 2.5 3y 0.95 0.84 0.86 1.06 1.5 0.72 1.9y' 1 1.5 2 2.5 3 3.5 4求质点在时刻1.8时的速度,并画出插值多项式的图像。
1)运用Hermite插值法画出图像,如图4-1,并求质点在时刻1.8时的速度。
>>clear>>clc>>X=[0.1 0.5 1 1.5 2 2.5 3;0.95 0.84 0.86 1.06 1.5 0.72 1.9;1 1.5 2 2.5 3 3.5 4];>> x=0.1:0.01:3;>> H=Hermite1(X,x);>> plot(x,H)>> hold on>> plot(X(1,:),X(2,:),'r*')>> H1_8=Hermite(X,1.8);>> plot(1.8,H1_8,'go')>> legend('插值图像','原始点','目标点');图4-1二、验证高次插值的Runge现象问题分析和算法设计(一)Language插值代码function [Ln] =Lagrange(X,x)%请输入2*n+1矩阵X,X中第一行每个元素都是插值节点,X中第二行每个元素都是插值节点对应的函数值;%第二章P24例一拉格朗日插值n=size(X,2);d=0;for m=1:1:nif x==X(1,m);d=m;breakendend运行结果和总结 运行结果 例:给定函数55,11)(2≤≤-+=x xx f ; (1) 验证表2-10的误差结果(高次插值的Runge 现象);(2) 以0.1为步长分别进行Language 插值、分段线性插值、分段三次Hermite插值,画出三种插值函数以及f(x)的图像,比较三种插值结果。
数值计算(数值分析)实验4-分段三次埃尔米特(hermite)插值【c程序实现+流程图】

实验四分段三次埃尔米特插值(一)实验目的掌握分段三次埃尔米特插值算法。
(二)实验项目内容1.写出计算步骤和流程图。
2.对每种算法分别用C或c#程序实现。
3.调试程序。
可用以下数据进行调试。
已知函数y=1/(1+x2)在区间[0,3]上取等距插值节点,求区间[0,3]上的分段三次埃尔米特插值函数,并利用它求出f(1.5)的近似值(0.3075)。
x0 1 2iy 1 0.5 0.2 iy 0 -0.5 -0.16 i(三)主要仪器设备微机(四)实验室名称公共计算机实验室(五)实验报告撰写实验四分段三次埃尔米特插值实验报告一、流程图二、 程序代码#include<stdio.h>#include<math.h>float f0(float x) N Y开始输入i x ,i y ,xy=0, j=0t=1i j ix x t t x x -=- i=0,…j-1,j+1,…n i y y ty =+j=n? 输出y结束j=j+1{return((x-1)*(x-1)*(2*x+1));}float f1(float x){return(x*x*(-2*x+3));}float g0(float x){return(x*(x-1)*(x-1));}float g1(float x){return(x*x*(x-1));}void main(){float x0,x1,x,y0,y1,yy0,yy1,h,p;printf("输入x0,x1,x,y0,y1和yy0,yy1的取值");scanf("%f%f%f%f%f%f%f",&x0,&x1,&x,&y0,&y1,&yy0,&yy1); h=x1-x0;p=y0*f0((x-x0)/h)+y1*f1((x-x0)/h)+h*yy0*g0((x-x0)/h)+h*yy1*g1((x-x0)/h);printf("%f\n",p);}三、运行结果【截图】。
数值分析实验,用程序实现Hermite插值法

《数值分析》实验报告实验序号:实验六 实验名称: Hermite 插值法1. 实验目的:学会Hermite 插值法,并应用该算法于实际问题.2. 实验内容:求一个函数ϕ(x )用来近似函数f (x ),用分段三次Hermit 插值的方法来求解近似函数ϕ(x )并画出近似函数图像及原函数图像。
设在区间[a,b]上,给定n+1个插值节点b x x x x a n =<<<<=...210和相应的函数值n y y y ,...,,10以及一阶导数值''1'0,...,,ny y y ,求一个插值函数)(x H ,满足以下条件: (1)),...,2,1,0()(,)(''n i y x H y x H i i i i === (2) )(x H 在每一个小区间[1,+j j x x ]上是三次多项式。
对于给定函数11-,2511)(2≤≤+=x x x f 。
在区间[]11-,上画出f (x )和分段三次Hermit 插值函数)(x H 的函数图像。
3. 实验分析:算法分析:1. 分段三次Hermit 插值的算法思想:分段三次Hermit 插值的做法是在每一个小区间上作三次Hermit 插值,因此在每一个插值节点上都需要构造两个插值基函数)(),(x H x h i i ,然后再作它们的线性组合。
分段三次Hermit 插值基函数如下:⎪⎩⎪⎨⎧≤≤----+=其它 0 ))(21()(1021010100x x x x x x x x x x x x h ⎪⎩⎪⎨⎧≤≤---=其它 0 ))(()(10210100x x x x x x x x x x H1,...,2,1 0 ))(21( ))(21()(1211112111-=⎪⎪⎪⎩⎪⎪⎪⎨⎧≤<----+≤≤----+=++++---n i x x x x x x x x x x x x x x x x x x x x x x x h i i i i i i i i i i-i i i i i i i 其它1,...,2,1 0 ))(( ))(()(12111211-=⎪⎪⎪⎩⎪⎪⎪⎨⎧≤<---≤≤---=+++--n i x x x x x x x x x x x x x x x x x x x H i i i i i i i i-i i i i i 其它⎪⎩⎪⎨⎧≤<----+=---其它 0 ))(21()(1-n 2111n n n n n n n n x x x x x x x x x x x x h ⎪⎩⎪⎨⎧≤<---=--其它 0 ))(()(1-n 211n n n n n n x x x x x x x x x x H 分段三次Hermit 插值函数是:∑=+=n i i i i i x H y x h y x H 0'))()(()( 4. 实验代码:// LDlg.cpp : implementation file//#include "stdafx.h"#include "L.h"#include "LDlg.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////////////////////////////////////////// CAboutDlg dialog used for App Aboutclass CAboutDlg : public CDialog{public:CAboutDlg();// Dialog Data//{{AFX_DATA(CAboutDlg)enum { IDD = IDD_ABOUTBOX };//}}AFX_DATA// ClassWizard generated virtual function overrides//{{AFX_VIRTUAL(CAboutDlg)protected:virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV support //}}AFX_VIRTUAL// Implementationprotected://{{AFX_MSG(CAboutDlg)//}}AFX_MSGDECLARE_MESSAGE_MAP()};CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD){//{{AFX_DATA_INIT(CAboutDlg)//}}AFX_DATA_INIT}void CAboutDlg::DoDataExchange(CDataExchange* pDX){CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(CAboutDlg)//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)//{{AFX_MSG_MAP(CAboutDlg)// No message handlers//}}AFX_MSG_MAPEND_MESSAGE_MAP()///////////////////////////////////////////////////////////////////////////// // CLDlg dialogCLDlg::CLDlg(CWnd* pParent /*=NULL*/): CDialog(CLDlg::IDD, pParent){//{{AFX_DATA_INIT(CLDlg)// NOTE: the ClassWizard will add member initialization here//}}AFX_DATA_INIT// Note that LoadIcon does not require a subsequent DestroyIcon in Win32m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);}void CLDlg::DoDataExchange(CDataExchange* pDX){CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(CLDlg)// NOTE: the ClassWizard will add DDX and DDV calls here//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(CLDlg, CDialog)//{{AFX_MSG_MAP(CLDlg)ON_WM_SYSCOMMAND()ON_WM_PAINT()ON_WM_QUERYDRAGICON()ON_BN_CLICKED(IDC_LARGRI, OnLargri)ON_BN_CLICKED(IDC_BUTTON2, OnButton2)ON_BN_CLICKED(IDC_HERMITE, OnHermite)//}}AFX_MSG_MAPEND_MESSAGE_MAP()///////////////////////////////////////////////////////////////////////////// // CLDlg message handlersBOOL CLDlg::OnInitDialog(){CDialog::OnInitDialog();// Add "About..." menu item to system menu.// IDM_ABOUTBOX must be in the system command range.ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);ASSERT(IDM_ABOUTBOX < 0xF000);CMenu* pSysMenu = GetSystemMenu(FALSE);if (pSysMenu != NULL){CString strAboutMenu;strAboutMenu.LoadString(IDS_ABOUTBOX);if (!strAboutMenu.IsEmpty()){pSysMenu->AppendMenu(MF_SEPARATOR);pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);}}// Set the icon for this dialog. The framework does this automatically // when the application's main window is not a dialogSetIcon(m_hIcon, TRUE); // Set big iconSetIcon(m_hIcon, FALSE); // Set small icon// TODO: Add extra initialization herereturn TRUE; // return TRUE unless you set the focus to a control}void CLDlg::OnSysCommand(UINT nID, LPARAM lParam){if ((nID & 0xFFF0) == IDM_ABOUTBOX){CAboutDlg dlgAbout;dlgAbout.DoModal();}else{CDialog::OnSysCommand(nID, lParam);}// If you add a minimize button to your dialog, you will need the code below // to draw the icon. For MFC applications using the document/view model, // this is automatically done for you by the framework.void CLDlg::OnPaint(){if (IsIconic()){CPaintDC dc(this); // device context for paintingSendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);// Center icon in client rectangleint cxIcon = GetSystemMetrics(SM_CXICON);int cyIcon = GetSystemMetrics(SM_CYICON);CRect rect;GetClientRect(&rect);int x = (rect.Width() - cxIcon + 1) / 2;int y = (rect.Height() - cyIcon + 1) / 2;// Draw the icondc.DrawIcon(x, y, m_hIcon);}else{CDialog::OnPaint();}}// The system calls this to obtain the cursor to display while the user drags // the minimized window.HCURSOR CLDlg::OnQueryDragIcon(){return (HCURSOR) m_hIcon;}void CLDlg::OnOK(){int x00=300,y00=350,i,j;double x;CDC *pDC=GetDC();pDC->SetMapMode(MM_LOMETRIC);pDC->SetViewportOrg(x00,y00);//画坐标轴与原函数for(i=-700; i<=700; i++){pDC->SetPixel(i,0,RGB(0,0,0));pDC->SetPixel(0,i,RGB(0,0,0));}for(x=-1; x<=1; x+=0.001){double j=400.0/(1+25*x*x);pDC->SetPixel(x*500,j,RGB(255,0,0));}pDC->TextOut(-30,-10,"0");pDC->TextOut(-30,430,"1");pDC->TextOut(490,-10,"1");pDC->TextOut(-490,-10,"-1");pDC->MoveTo(-10,680); //x箭头pDC->LineTo(0,700);pDC->MoveTo(0,700);pDC->LineTo(10,680);pDC->MoveTo(680,10); //y箭头pDC->LineTo(700,0);pDC->MoveTo(700,0);pDC->LineTo(680,-10);pDC->TextOut(-30,700,"y");pDC->TextOut(700,-10,"x");}void CLDlg::OnLargri(){int x00=300,y00=350,i,j;double x;CDC *pDC=GetDC();pDC->SetMapMode(MM_LOMETRIC);pDC->SetViewportOrg(x00,y00);//画坐标轴for(i=-700; i<=700; i++){pDC->SetPixel(i,0,RGB(0,0,0));pDC->SetPixel(0,i,RGB(0,0,0));}double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1}; pDC->TextOut(-30,-10,"0");pDC->TextOut(-30,430,"1");pDC->TextOut(490,-10,"1");pDC->TextOut(-490,-10,"-1");pDC->MoveTo(-10,680); //x箭头pDC->LineTo(0,700);pDC->MoveTo(0,700);pDC->LineTo(10,680);pDC->MoveTo(680,10); //y箭头pDC->LineTo(700,0);pDC->MoveTo(700,0);pDC->LineTo(680,-10);pDC->TextOut(-30,700,"y");pDC->TextOut(700,-10,"x");// 拉格朗日差值的函数double yy[12],lx[12],ly[12];double l_fenzi[12],l_fenmu[12];double l_x,l_y;for(i=0; i<=10; i++){yy[i]=1.0/(1+25*yx[i]*yx[i]);for(i=0; i<=10; i++){l_fenmu[i]=1.0;for(j=0; j<=10; j++){if(i!=j)l_fenmu[i]=l_fenmu[i]*(yx[i]-yx[j]);}}double qq,pp;for(qq=-1; qq<=1; qq+=0.0003){for(i=0; i<=10; i++){l_fenzi[i]=1.0;for(j=0; j<=10; j++){if(i!=j)l_fenzi[i]=l_fenzi[i]*(qq-yx[j]);}}pp=0;for(i=0; i<=11; i++){pp=pp+1.0/(1+25*yx[i]*yx[i])*l_fenzi[i]/l_fenmu[i];}pDC->SetPixel(qq*500,pp*390+5,RGB(132,112,225));}void CLDlg::OnButton2(){int x00=300,y00=350,i,j;double x;CDC *pDC=GetDC();pDC->SetMapMode(MM_LOMETRIC);pDC->SetViewportOrg(x00,y00);//画坐标轴与原函数for(i=-700; i<=700; i++){pDC->SetPixel(i,0,RGB(0,0,0));pDC->SetPixel(0,i,RGB(0,0,0));}double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1}; double yy[14];for(i=0; i<=10; i++){yy[i]=1.0/(1+25*yx[i]*yx[i]);}pDC->TextOut(-30,-10,"0");pDC->TextOut(-30,430,"1");pDC->TextOut(490,-10,"1");pDC->TextOut(-490,-10,"-1");pDC->MoveTo(-10,680); //x箭头pDC->LineTo(0,700);pDC->MoveTo(0,700);pDC->LineTo(10,680);pDC->MoveTo(680,10); //y箭头pDC->LineTo(700,0);pDC->MoveTo(700,0);pDC->LineTo(680,-10);pDC->TextOut(-30,700,"y");pDC->TextOut(700,-10,"x");// 线性分段差值的图像CPen pen;CPen*oldpen;pen.CreatePen(PS_SOLID,5,RGB(0,0,0));oldpen=pDC->SelectObject(&pen);for(i=0; i<10; i++){pDC->MoveTo(yx[i]*480,yy[i]*400);pDC->LineTo(yx[i+1]*480,yy[i+1]*400); }}void CLDlg::OnHermite(){int x00=300,y00=350,i,j;double x;CDC *pDC=GetDC();pDC->SetMapMode(MM_LOMETRIC);pDC->SetViewportOrg(x00,y00);//画坐标轴与原函数for(i=-700; i<=700; i++){pDC->SetPixel(i,0,RGB(0,0,0));pDC->SetPixel(0,i,RGB(0,0,0));}double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1};double yy[12];for(i=0; i<=10; i++){yy[i]=1.0/(1+25*yx[i]*yx[i]);}pDC->TextOut(-30,-10,"0");pDC->TextOut(-30,430,"1");pDC->TextOut(490,-10,"1");pDC->TextOut(-490,-10,"-1");pDC->MoveTo(-10,680); //x箭头pDC->LineTo(0,700);pDC->MoveTo(0,700);pDC->LineTo(10,680);pDC->MoveTo(680,10); //y箭头pDC->LineTo(700,0);pDC->MoveTo(700,0);pDC->LineTo(680,-10);pDC->TextOut(-30,700,"y");pDC->TextOut(700,-10,"x");//分段三次Hermite差值的函数double x0,x1,yd1,yd0,y1,y0;for(i=0; i<10; i++){x0=yx[i],x1=yx[i+1];y0=1.0/(1+25*x0*x0);y1=1.0/(1+25*x1*x1);yd0=-(50*x0)*1.0/((1+25*x0*x0)*(1+25*x0*x0));yd1=-(50*x1)*1.0/((1+25*x1*x1)*(1+25*x1*x1));for(double qq=x0; qq<x1; qq+=0.005){double pp= y0*(1+2*(qq-x0)/(x1-x0)) * (qq-x1)/(x0-x1) * (qq-x1)/(x0-x1)+y1*(1+2*(qq-x1)/(x0-x1)) * (qq-x0)/(x1-x0) * (qq-x0)/(x1-x0)+yd0*(qq-x0) * (qq-x1)/(x0-x1) * (qq-x1)/(x0-x1)+yd1*(qq-x1) * (qq-x0)/(x1-x0) * (qq-x0)/(x1-x0);pDC->SetPixel(qq*500,pp*400,RGB(225,185,15));}}}5.实验截图6. 实验结果分析:通过本次实验我对分段三次Hermit插值有了更深刻更全面的掌握,它在给定了节点处的函数值和导数值以后,构造了一个整体上具有一阶连续微商的插值函数。
分段三次埃尔米插值

实习报告题目分段三次埃尔米插值专业信息与计算科学班级计算092学号3090811050学生杜青2011 年一、目的意义:分段线性插值多项式在插值区间上只能保证连续性,而不光滑。
要想得到在差值区间上光滑的分段差值多项式,可采用分段三次埃尔米差值多项式。
二、 数学公式22111332211122[2()]()[2()]()()()()()()i i i i i i i i i i ii i i i i i i ih x x x x h x x x x H x y y h h x x x x x x x x y y h h ------+-----=+----''++1012121(),[,](),[,]()(),[,]n n n H x x x x H x x x x H x H x x x x -∈⎧⎪∈⎪=⎨⎪⎪∈⎩ 三、 算法#include <stdio.h>void main(){ int i;float s,X,m1,m2,m3,m4,h;float x[20],y[20];printf("请输入要插入的值:x=");scanf("%f",&X);for(i=0;i<4;i++){printf("x[%d]=",i);scanf("%f",&x[i]);printf("y[%d]=",i);scanf("%f",&y[i]);}h=x[1]-x[0];m1=((x[1]-x[0])+2*(X-x[0]))*(X-x[1])*(X-x[1])*y[0]/h/h/h;m2=((x[1]-x[0])-2*(X-x[1]))*(X-x[0])*(X-x[0])*y[1]/h/h/h;m3=(X-x[0])*(X-x[1])*(X-x[1])*y[2]/h/h;m4=(X-x[0])*(X-x[0])*(X-x[1])*y[3]/h/h;s=m1+m2+m3+m4;printf("%f",s);}四、 流程图五、算法实例已知(1,2),(2,3),(1,1),(2,-1),插入的值为1.5,求所得的结果。
数值分析(13)Hermite插值

a 2li' ( xi ) 解出 ' b 1 2 x l i i ( xi )
hi ( x ) (1 2( x xi )li' ( xi )) l i2 ( xi ) ( i 0,1, 2, n) n n x xj 1 ' 其中 li ( x ) ( ), l i ( xi ) ( ) xi x j xi x j j0 j0 所以
2n
共2n 2个方程,可求出2n 2个系数a0 , a1 ,..., a2n , a2n1 .
数值分析
数值分析
Hermite插值多项式的构造
( 2) Lagrange型插值基函数法 设Hermite插值多项式为 H 2 n1 ( x ) hi ( x ) yi hi ( x ) y 'i
数值分析
数值分析
由条件(2)可列出方程组 2 h ( x ) ( cx d ) l i i i i ( xi ) 0 2 ' h ' ( x ) cl ( x ) 2( cx d ) l ( x ) l i i i i i i i i ( xi ) 1
li ( xi ) 1, cxi d 0, c 1 解出 d xi 于是求出
i 0
数值分析
n
2
F ( t )关于t 有n 2个零点:x0,x1, ,xn,x 。 但F ' ( t )关于t 有2n 2个零点,由Rolle(罗尔)定理 必存在点 (a , b),使 F
(2 n 2)
( ) f
(2 n 2)
( ) 0 K ( x )(2n 2)! 0
实习:Matlab作业hermite插值

题目:利用Matlab实现数据的Hermite插值和分段三次Hermite插值小组成员:王晓波(38)蔡明宇(20)一、程序实现意义:一般的,从各种试验得来的数据总有一定的数量,而利用插值技术能够从有限的数据中获取整体的状态。
而Hermite插值不仅保证了插值函数与原函数在给定数据点处得拟合,同时保证了在相应点处导数的相同,从而在很大程度上保证了曲线的“光滑性”。
因此,通过Matlab实现Hermite插值具有很普遍的意义。
二、实现过程:1、Hermite插值由于并不是所有的Matlab版本都提供现有的Hermite插值函数包,故我们首先编写了实现给定五个观测点的Hermite插值的M程序,代码如下:function [f,f0] = Hermite1(x,y,y_1)syms t;f = ;!if(length(x) == length(y))if(length(y) == length(y_1))n = length(x);elsedisp('y和y的导数的维数不相等');return;endelsedisp('x和y的维数不相等! ');return;end*for i=1:nh = ;a = ;for j=1:nif( j ~= i)h = h*(t-x(j))^2/((x(i)-x(j))^2);a = a + 1/(x(i)-x(j));endendf = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));<endf0 = subs(f,'t');其中x为给定点横坐标数组,y为给定点纵坐标数组,y_1为原函数在给定点处的导数数组。
测试证明该程序可以实现,例如输入如下数组:x=1::;y_1=[ ];y=[1 ];>> [f,f0] = Hermite1(x,y,y_1);运行结果如下:f =$(390625*((3972231*t)/35 - 28321/0000)*(t - 1)^2*(t - 7/5)^2*(t - 8/5)^2*(t - 9/5)^2)/36 - (390625*(t - 1)^2*(t - 6/5)^2*(t - 7/5)^2*(t - 9/5)^2*((28557*t)/28 - 23/2000))/36 + (390625*((64*t)/3 - 61/3)*(t - 6/5)^2*(t - 7/5)^2*(t - 8/5)^2*(t - 9/5)^2)/576 + (390625*((763*t)/1984 + 043/6240000)*(t - 1)^2*(t - 6/5)^2*(t - 8/5)^2*(t - 9/5)^2)/16 - (390625*((77623*t)/28 - 931/60000)*(t - 1)^2*(t - 6/5)^2*(t - 7/5)^2*(t - 8/5)^2)/576f0 =.利用matlab绘制图像:2、程序的窗口化:利用Matlab提供的GUIDE工具以及callback函数实现相应函数的窗口化,GUI代码如下:function varargout = untitled(varargin)?% UNTITLED M-file for% UNTITLED, by itself, creates a new UNTITLED or raises the existing% singleton*.%% H = UNTITLED returns the handle to a new UNTITLED or the handle to% the existing singleton*.%% UNTITLED('CALLBACK',hObject,eventData,handles,...) calls the local% function named CALLBACK in with the given input arguments.%% UNTITLED('Property','Value',...) creates a new UNTITLED or raises the,% existing singleton*. Starting from the left, property value pairs are% applied to the GUI before untitled_OpeningFcn gets called. An% unrecognized property name or invalid value makes property application% stop. All inputs are passed to untitled_OpeningFcn via varargin.%% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one% instance to run (singleton)".%% See also: GUIDE, GUIDATA, GUIHANDLES% Edit the above text to modify the response to help untitled%% Last Modified by GUIDE 15-Sep-2011 22:24:48% Begin initialization code - DO NOT EDITgui_Singleton = 1;gui_State = struct('gui_Name', mfilename, ...'gui_Singleton', gui_Singleton, ...'gui_OpeningFcn', @untitled_OpeningFcn, ...'gui_OutputFcn', @untitled_OutputFcn, ...'gui_LayoutFcn', [] , ...'gui_Callback', []);】if nargin && ischar(varargin{1})= str2func(varargin{1});endif nargout[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); elsegui_mainfcn(gui_State, varargin{:});end% End initialization code - DO NOT EDIT<% --- Executes just before untitled is made visible.function untitled_OpeningFcn(hObject, eventdata, handles, varargin) % This function has no output args, see OutputFcn.% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)% varargin command line arguments to untitled (see VARARGIN)% Choose default command line output for untitled= hObject;>% Update handles structureguidata(hObject, handles);% UIWAIT makes untitled wait for user response (see UIRESUME)% uiwait;% --- Outputs from this function are returned to the command line.function varargout = untitled_OutputFcn(hObject, eventdata, handles)% varargout cell array for returning output args (see VARARGOUT);…% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Get default command line output from handles structurevarargout{1} = ;function edit1_Callback(hObject, eventdata, handles)% hObject handle to edit1 (see GCBO)…% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of edit1 as text% str2double(get(hObject,'String')) returns contents of edit1 as a double guidata(hObject, handles);% --- Executes during object creation, after setting all properties.function edit1_CreateFcn(hObject, eventdata, handles)% hObject handle to edit1 (see GCBO)(% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end¥function edit2_Callback(hObject, eventdata, handles)% hObject handle to edit2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of edit2 as text% str2double(get(hObject,'String')) returns contents of edit2 as a double guidata(hObject, handles);% --- Executes during object creation, after setting all properties.、function edit2_CreateFcn(hObject, eventdata, handles)% hObject handle to edit2 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end—function edit3_Callback(hObject, eventdata, handles)% hObject handle to edit3 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of edit3 as text% str2double(get(hObject,'String')) returns contents of edit3 as a double guidata(hObject, handles);·% --- Executes during object creation, after setting all properties.function edit3_CreateFcn(hObject, eventdata, handles)% hObject handle to edit3 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');·endfunction edit4_Callback(hObject, eventdata, handles)% hObject handle to edit4 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of edit4 as text% str2double(get(hObject,'String')) returns contents of edit4 as a double '% --- Executes during object creation, after setting all properties.function edit4_CreateFcn(hObject, eventdata, handles)% hObject handle to edit4 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'));set(hObject,'BackgroundColor','white');end% --- Executes on button press in pushbutton1.function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)x=str2num(get,'string'));y=str2num(get,'string'));<y_1=str2num(get,'string'));x0=str2num(get,'string'));syms t;f = ;if(length(x) == length(y))if(length(y) == length(y_1))n = length(x);elsedisp('yºÍyµÄµ¼ÊýµÄάÊý²»ÏàµÈ');return;end—elsedisp('xºÍyµÄάÊý²»ÏàµÈ£¡ ');return;endfor i=1:nh = ;a = ;for j=1:nif( j ~= i)h = h*(t-x(j))^2/((x(i)-x(j))^2);a = a + 1/(x(i)-x(j));、endendf = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));endf0 = subs(f,'t',x0);plot,x,y,'*');^function edit5_Callback(hObject, eventdata, handles)% hObject handle to edit5 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)% Hints: get(hObject,'String') returns contents of edit5 as text% str2double(get(hObject,'String')) returns contents of edit5 as a doubleguidata(hObject, handles);{% --- Executes during object creation, after setting all properties.function edit5_CreateFcn(hObject, eventdata, handles)% hObject handle to edit5 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles empty - handles not created until after all CreateFcns called% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc && isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))~set(hObject,'BackgroundColor','white');end程序运行结果:其中左上方纵列的三个对话框从上到下分别输入给定点的横坐标x,纵坐标y以及导数值y_1,右侧空白框输入维数,下方坐标图显示插值函数图像,例如仍插入上面所给定的点列,得出结果:从图上看拟合程度还是比较不错的。
分段三次Hermite插值计算GNSS系统卫星仰角

分段三次Hermite插值计算GNSS系统卫星仰⾓龙源期刊⽹ /doc/847794544.html分段三次Hermite插值计算GNSS系统卫星仰⾓作者:李路程来源:《科技视界》2016年第15期【摘要】⽂中采⽤300s时间间隔的分段三次Hermite插值多项式计算GNSS系统卫星仰⾓,通过与理论仰⾓的⽐较和精度统计,验证⽂中⽅法的正确性。
【关键词】GNSS系统;卫星仰⾓;分段三次Hermite插值0 引⾔在全球卫星定位技术⽅兴未艾之际,随着越来越多新技术、新产品的研制开发,在产品定型之前,需要尽可能的发现潜在的问题,确保技术产品的可靠性。
试验条件下事后数据处理⼯作的重要性逐渐受到⼈们的关注和认可。
数据分析⼯作中经常需要筛选出⾼度⾓在条件区间的数据信息,以衡量接收机数据质量的指标,如伪距精度、载波精度、完整率等等。
本⽂的⼯作定位于根据接收机记录的⾼频观测数据和星历⽂件,计算卫星的仰⾓,GNSS系统选择为GPS、GLONASS、Galileo三个系统。
3 算例计算试验过程中接收机接收模拟源数字信号,记录星历⽂件、1Hz的观测数据,数据时长17000s。
算例中将⽂中⽅法根据记录的观测数据、星历⽂件解算的卫星仰⾓,和模拟源仿真数据的卫星仰⾓作⽐对,验证⽂中⽅法在GNSS系统卫星仰⾓计算的正确性。
下⾯依次给出了GPS、GLONASS和Galileo三个系统的卫星仰⾓对⽐图。
为便于作图⽐较,先将⽂中⽅法计算的卫星S的仰⾓在数值上增加1.0,再与模拟源记录的理论仰⾓信息做⽐对图。
图表中MNY_S*表⽰模拟源的卫星*的仰⾓,S*+1.0表⽰⽂中⽅法计算的卫星*的仰⾓数值增加1.0。
GPS系统选取的卫星为03、16、32,三颗卫星的仰⾓信息如图1所⽰。
图1表明⽂中⽅法计算的GPS仰⾓与模拟源的变化过程相符,能真实地反映出卫星起落现象,同时证实接收机在卫星仰⾓低于10度时,依然能够完好地进⾏信号的跟踪捕获。
GLONASS系统选取的卫星为02、14、24,三颗卫星的仰⾓信息如图2所⽰。
分段三次Hermite插值(可编辑修改word版)

摘要用函数来表示变量间的数量关系广泛应用于各学科领域,但是在实际问题中,往往是通过实验、观测以及计算等方法,得到的是函数在一些点上的函数值。
如何通过这些离散数据找到函数的一个满足精度要求且便于使用的近似表达式,是经常遇到的问题。
对于这类问题我们解决的方法为插值法,而最常用也最简单的插值方法就是多项式插值。
当然用插值法得到的近似表达式必须满足插值条件即假设给定了 n+1 个点的自变量的值以及函数值,近似函数必须要过这n+1 个点。
多项式插值,从几何角度看,就是寻求 n 次代数曲线y=Pn (x)通过 n+1 个点作为 f(x)的近似。
但是随着插值节点个数的增加,高次插值多项式的近似效果并不理想。
根据大量实验得出,在进行高次多项式插值时,会出现龙格现象。
因此,为了解决这样的一个问题,我们可以通过缩小插值区间的办法达到减小误差的目的。
但是当在每个小区间上用一次函数进行插值时,有很好的收敛性但是光滑度不够,因此本实验将用三次 Hermite 进行插值,做具体的讨论和学习。
关键词:龙格现象分段差值三次Hermite 进行插值i i ‒ 1 i ‒ 1 i iH 2(x ),x ∈ [x 1,x 2] ⋯⋯⋯⋯⋯⋯⋯ H n (x ),x ∈ [x n ‒ 1,x n ]1、实验目的1) 通过对分段三次 Hermite 插值算法程序的编写,提高自己编写程序的能力2) 体会分段三次 Hermite 插值比分段线性插值优越在哪里 3) 用实验报告的形式展现,提高自己在写论文方面的能力2、算法流程分段线性插值多项式S (x )在插值区间[a ,b ]上只能保证连续性,而不光滑。
要想得到在插值区间上光滑的分段线性插值多项式,可采用分段埃尔米特(Hermite )插值,这里我们考虑在整个[a ,b ]上用分段三次埃尔米特插值多项式来逼近f (x )。
一般的将带有导数的插值多项式称为 Hermite 插值多项式。
如果已知函数y = f (x )在节点a = x 0 < x 1 < ⋯ < x n = b 处的函数的值 和导数值:y i = f (x i ),y ' = f '(x i ),i = 0,1,2,⋯,n则在小区间[x i ‒ 1,x i ]上有四个插值条件:y i ‒ 1 = f (x i ‒ 1),y i = f (x i )y ' = f '(x ),y ' = f '(x ) 故能构造一个三次多项式H i (x ),并称为三次 Hermite 插值多项式。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验四分段三次埃尔米特插值
(一)实验目的
掌握分段三次埃尔米特插值算法。
(二)实验项目内容
1.写出计算步骤和流程图。
2.对每种算法分别用C或c#程序实现。
3.调试程序。
可用以下数据进行调试。
已知函数y=1/(1+x2)在区间[0,3]上取等距插值节点,求区间[0,3]上的分段三次埃尔米特插值函数,并利用它求出f(1.5)的近似值(0.3075)。
x0 1 2
i
y 1 0.5 0.2 i
y 0 -0.5 -0.16 i
(三)主要仪器设备
微机
(四)实验室名称
公共计算机实验室
(五)实验报告撰写
实验四分段三次埃尔米特插值
实验报告
一、流程图
二、 程序代码
#include<stdio.h>
#include<math.h>
float f0(float x) N Y
开始
输入i x ,i y ,x
y=0, j=0
t=1
i j i
x x t t x x -=- i=0,…j-1,j+1,…n i y y ty =+
j=n? 输出y
结束
j=j+1
{
return((x-1)*(x-1)*(2*x+1));
}
float f1(float x)
{
return(x*x*(-2*x+3));
}
float g0(float x)
{
return(x*(x-1)*(x-1));
}
float g1(float x)
{
return(x*x*(x-1));
}
void main()
{float x0,x1,x,y0,y1,yy0,yy1,h,p;
printf("输入x0,x1,x,y0,y1和yy0,yy1的取值");
scanf("%f%f%f%f%f%f%f",&x0,&x1,&x,&y0,&y1,&yy0,&yy1); h=x1-x0;
p=y0*f0((x-x0)/h)+y1*f1((x-x0)/h)+h*yy0*g0((x-x0)/h)+h*yy1*g1((x-
x0)/h);
printf("%f\n",p);
}
三、运行结果【截图】。