hermite插值以及两种MATLAB程序
课程设计---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 插值问题的一般性公式。
计算方法-插值方法实验
实验一插值方法一. 实验目的（1）熟悉数值插值方法的基本思想，解决某些实际插值问题，加深对数值插值方法的理解。
（2）熟悉Matlab 编程环境，利用Matlab 实现具体的插值算法，并进行可视化显示。
二. 实验要求用Matlab 软件实现Lagrange 插值、分段线性插值、三次Hermite 插值、Aitken 逐步插值算法，并用实例在计算机上计算和作图。
三. 实验内容1. 实验题目 （1）已知概率积分dxe y xx ⎰-=22π的数据表构造适合该数据表的一次、二次和三次Lagrange 插值公式，输出公式及其图形，并计算x =0.472时的积分值。
答：①一次插值公式：输入下面内容就可以得到一次插值结果 >> X=[0.47,0.48];Y=[0.4937452,0.5027498]; >> x=0.472;>> (x-X(2))/(X(1)-X(2))*Y(1)+(x-X(1))/(X(2)-X(1))*Y(2)ans =0.495546120000000>>②两次插值公式为：输入下面内容就可以得到两次插值结果>> X=[0.46,0.47,0.48];Y=[0.4846555,0.4937452,0.5027498]; >> x=0.472;>>(x-X(2))*(x-X(3))/((X(1)-X(2))*(X(1)-X(3)))*Y(1)+(x-X(1))*(x-X(3))/((X(2)-X(1))*(X(2)-X(3)))*Y(2)+(x-X(2))*(x-X(1))/((X(3)-X(2))*(X(3)-X(1)))*Y(3)i 0123x 0.46 047 0.48 0.49 y0.4846555 0.4937452 0.5027498 0.5116683ans =0.495552928000000>>③三次插值公式为：输入下面内容就可以得到三次插值结果>> X=[0.46,0.47,0.48,0.49];Y=[0.4846555,0.4937452,0.5027498,0.5116683];>> x=0.472;>>(x-X(2))*(x-X(3))*(x-X(4))/((X(1)-X(4))*(X(1)-X(2))*(X(1)-X(3)))*Y(1)+(x-X(4))*( x-X(1))*(x-X(3))/((X(2)-X(4))*(X(2)-X(1))*(X(2)-X(3)))*Y(2)+(x-X(4))*(x-X(2))*( x-X(1))/((X(3)-X(4))*(X(3)-X(2))*(X(3)-X(1)))*Y(3)+(x-X(3))*(x-X(2))*(x-X(1))/(( X(4)-X(1))*(X(4)-X(2))*(X(4)-X(3)))*Y(4)ans =0.495552960000000输入下面内容，绘出三点插值的图：>> X=[0.46,0.47,0.48,0.49];Y=[0.4846555,0.4937452,0.5027498,0.5116683];>> x=linspace(0.46,0.49);>>y=(x-X(2)).*(x-X(3)).*(x-X(4))/((X(1)-X(4))*(X(1)-X(2))*(X(1)-X(3)))*Y(1)+(x-X(4) ).*(x-X(1)).*(x-X(3))/((X(2)-X(4))*(X(2)-X(1))*(X(2)-X(3)))*Y(2)+(x-X(4)).*(x-X(2) ).*(x-X(1))/((X(3)-X(4))*(X(3)-X(2))*(X(3)-X(1)))*Y(3)+(x-X(3)).*(x-X(2)).*(x-X(1) )/((X(4)-X(1))*(X(4)-X(2))*(X(4)-X(3)))*Y(4);>>plot(x,y)（注意上面的“.*”不能用“*”替代）；（2）将区间[-5,5]分为10等份，求作211)(x x f +=的分段线性插值函数，输出函数表达式及其图形，并计算x =3.3152时的函数值。
埃尔米特(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; }。
5.2Hermite插值
六、分段Hermite插值 （与学生一起看书学习） 6.6 分段Hermite插值（任玉杰） 1、分段Hermite插值函数（408页） （1）定义6.4（分段三次Hermite插值函数）
分段三次Hermite插值函数与一般的Hermite插值的区别：
设在节点 a x0 x1 xn b 上,
r ( xi ) 0, i 0,1,2,, n r ( xik ) 0, k 0,1,2,, m
r ( x) 的零点个数m+n+2个.
r ( x) 0
四、Hermite插值多项式的构造 设在节点 a x0 x1 xn b 上,
yi f xi , yi f xi i 0,1, n
0 ( x) y1 1 ( x) H3 ( x) y00 ( x) y11 ( x) y0
H 2 n1 x y j j x yj j x .
j 0 n
n 1 2 j x 1 2 x x j l j x. k 0 x j xk k j
提示：令
j ( x) (ax b)l 2 j ( x), j ( x) (cx d )l 2 j ( x),
例：求满足条件：
f x0 , y1 f x1 , y1 f x1 y0 f x0 , y0
的三次Hermite插值多项式 H 3 ( x) 。
要求插值多项式 H x ,满足条件
① H xi yi , H xi yi i 0,1, n .
1 H x C [a, b] ②
③在每个子区间 [ xi 1 , xi ] 上为三次多项式。 （2）分段三次Hermite插值多项式的形式
(完整版)Matlab学习系列13.数据插值与拟合
13. 数据插值与拟合实际中，通常需要处理实验或测量得到的离散数据（点）。
插值与拟合方法就是要通过离散数据去确定一个近似函数（曲线或曲面），使其与已知数据有较高的拟合精度。
1.如果要求近似函数经过所已知的所有数据点，此时称为插值问题（不需要函数表达式）。
2.如果不要求近似函数经过所有数据点，而是要求它能较好地反映数据变化规律，称为数据拟合（必须有函数表达式）。
插值与拟合都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数。
区别是：【插值】不一定得到近似函数的表达形式，仅通过插值方法找到未知点对应的值。
【拟合】要求得到一个具体的近似函数的表达式。
因此，当数据量不够，但已知已有数据可信，需要补充数据，此时用【插值】。
当数据基本够用，需要寻找因果变量之间的数量关系（推断出表达式），进而对未知的情形作预测，此时用【拟合】。
一、数据插值根据选用不同类型的插值函数，逼近的效果就不同，一般有：（1）拉格朗日插值（lagrange插值）（2）分段线性插值（3）Hermite（4）三次样条插值Matlab 插值函数实现：（1）interp1( ) 一维插值（2）intep2( ) 二维插值（3）interp3( ) 三维插值（4）intern( ) n维插值1.一维插值（自变量是1维数据）语法：yi = interp1(x0, y0, xi, ‘method’)其中，x0, y0为原离散数据（x0为自变量，y0为因变量）；xi为需要插值的节点，method为插值方法。
注：（1）要求x0是单调的，xi不超过x0的范围；（2）插值方法有‘nearest’——最邻近插值；‘linear’——线性插值；‘spline’——三次样条插值；‘cubic’——三次插值；默认为分段线性插值。
例1 从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24．试估计每隔1/10小时的温度值。
埃尔米特(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; }。
matlab实现newton差值和hermite差值
（一）实验目的掌握并能够利用newton差值和hermite差值方法解决问题。
（二）问题描述问题四插值。
上述函数的导数为采用三种方法中最好的方法计算这一积分(1)利用数值积分的方法给出在(可以直接计算精确值的，用精确值)，用Newton插值方法得到5个椭圆的周长(2)利用数值积分的方法给出在(可以直接计算精确值的，用精确值)，用Hermite插值方法得到5个椭圆的周长(3) 选做题：利用以及导数更多的值来进行插值，插值误差会有什么变化？(4)选做题：采用其它的插值方法改进插值的效果。
(三)算法介绍a确定，对于给定的b值都对应着一个椭圆，在本问题中用newton插值法和hermite得到的多项式代替椭圆周长公式中的进行积分，首先画出图像，选择初始点。
图像的实现代码见picture1.m。
newton插值法迭代公式：；Hermite法迭代公式：。
（四）程序建立picture.m文件画出和其导数图像。
（注：此图像为b=0.5时）x=0:0.1:2;y=sqrt(1+(0.5^2-1).*cos(x).^2);yyy=.750./(1-.75.*cos(x).^2).^(1/2).*cos(x).*sin(x);plot(x,y,'r');hold on;plot(x,yyy);hold off;legend('sqrt(1+(0.5^2-1).*cos(x).^2)','.750./(1-.75.*cos(x).^2).^(1/2).*cos(x).*sin(x)');所画图像为：我们选取0,0.3,0.6,0.9,1.2,1.5为初始点。
问题四（1）建立newtondedai1.m文件。
function z=newtondedai1(f,n)syms xia=zeros(n,n);x=[0 0.3 0.6 0.9 1.2 1.5];y=feval(f,x);a(:,1)=y;for i=2:nfor j=2:ia(i,j)=(a(i,j-1)-a(i-1,j-1))/(x(1,i)-x(1,i-j+1)); endendt=xi-x(1,1);p=a(1,1);for i=2:np=p+a(i,i)*t;t=t*(xi-x(1,i));endp=collect(vpa(p))问题四（2）建立hermite3.m文件。
第2章(2)Hermite插值
根据 j (x)及 j (x) 的一般表达式（4.4）及（4.5）， 可得到
x xk x xk 1 k ( x ) 1 2 x x , xk 1 xk k k 1 2 x xk 1 x xk k 1 ( x) 1 2 x x . xk xk 1 k 1 k
2.4
Hermite插值
【应用中进一步的插值问题】
（1）要求在节点上函数值相等；
（2）对应的导数值也相等，甚至要求高阶导数也相等. 满足这种要求的插值多项式就是Hermite插值多项式. 下面只讨论函数值与导数值个数相等的情况.
1
2.4.1
重节点均差与Taylor插值
定理3 设 f C n [a, b], x0 , x1 ,, xn 为 [a, b] 上的相异节点，
12
两端取对数再求导，得
l ( x j ) j
k 0 k j
n
1 , x j xk
于是
j ( x) (1 2( x x j )
k 0 k j n
1 )l 2 ( x). j x j xk
(4.4)
同理，可得
j ( x) ( x x j )l 2 ( x). j
2
2
同理可得
7
再求 0( x) ，可设
x x1 0 ( x) a( x x0 ) x x , 1 0
2
显然
直接由
0 ( x1 ) 0 ' ( x1 ) 0.
0 ' ( x0 ) a 1 解得
x x1 0 ( x) ( x x0 ) x x . 1 0
第三章 插值法 Hermite插值
j
(
x
)
(1
c(
x
x
j
))
(((xxx xx000))222((xx xx11))222 ((xxjj xx00))22((xxjj xx11))22
(((xxxxxjjj11))222((xx xxjj11))22 ((xx ((xxjj xxjj11))22((xxjj xxjj11))22 ((
P3( x) f ( x0 ) f x0, x1x x0 f x0, x1, x2 ( x x0 )x x1
f [ x0 , x1, x1, x2 ,]( x x0 )( x x1 )( x x2 ) --- 带重节点的牛顿插值多项式
插值余项（误差估计）:
条件 f ( x) C 3[a, b]，f（4）( x)存在。
i 0,1, , n, xi互异)
（b）H2n1( x) 为Hermite插值多项式，
则
R2n1( x) f ( x) H 2n1( x)
f ( (2n2) )
(2n 2)!
(x
x0 )2( x
x1 )2
( x xn )2
f (2n (2n
2) ( )
2)!
2 n1
(
x
),
(a, b)且与 x 有关。
结论
R( x)
f ( x) P3 ( x)
f
(4) (
4!
)
(
x
x0
)( x
x1 )2 ( x x0 , x2
x2 )。 且依赖于x
证明
（1）当 x xi (i 0,1,2)时, R( xi ) 0 右端 0; （2）当 x xi (i 0,1,2)时,
MATLAB应用多项式插值
多项式为 y 0.2015x3 1.4385x2 2.7477x 5.4370 Polyfit的第三个参数是多项式的阶数。
要求导数值也相等，甚至要求高阶导数值也相等，满足 这一要求的插值多项式就是Hermite插值多项式。下面 只讨论函数值与一阶导数值个数相等且已知的情况。
已知n个插值点 x1, x2, , xn 及对应的函数值 y1, y2, , yn 和一阶导数值 y1', y2', , y'n 。则对插值区间 内任意x的函数值y的Hermite插值公式：
n
y(x) hi[(xi x)(2ai yi yi' ) yi ] i 1
其中
hi
n ( x xj )2; j1 xi x j
ji
ai
n j 1
xi
1 xj
ji
• MATLAB实现
% hermite.m
function y=hermite(x0,y0,y1,x)
n
y(x) hi[(xi x)(2ai yi yi' ) yi ]
• polyval: 可用命令polyval计算多项式的值。 例： y 3x4 7x3 2x2 x 1 计算y(2.5)
>> c=[3,-7,2,1,1]； xi=2.5; yi=polyval(c,xi) yi =
23.8125 如果xi是含有多个横坐标值的数组，则yi也
为与xi长度相同的向量。 >> c=[3,-7,2,1,1]; xi=[2.5,3]; >> yi=polyval(c,xi) yi =
数值分析实验报告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)的图像，比较三种插值结果。
分段Hermite插值
6.6 分段埃尔米特插值及其MATLAB 程序6.6.2 分段埃尔米特插值的MATLAB 程序调用格式一：YI=interp1(X,Y,XI,'pchip')调用格式二：YI=interp1(X,XI,'pchip')例6.6.5 试用MA TLAB 程序计算例6.6.1中在各小区间中点处分段三次埃尔米特插值)(2/1+i n x H 及其相对误差.解 在MATLAB 工作窗口输入程序>> h=0.2;x0=-1:h:1;y0=1./(1+25.*x0.^2); xi=-0.9:h:0.9; fi=1./(1+25.*xi.^2); yi=interp1(x0,y0,xi,'pchip'); Ri=abs((fi-yi)./fi); xi,fi,yi,Ri,i=[xi',fi',yi',Ri'] 运行后屏幕显示各小区间中点x i 处的函数值f i ，插值s i ，相对误差值R i (略).6.6.3 作有关分段埃尔米特插值图形的MATLAB 程序作有关分段埃尔米特插值图形的MATLAB 主程序function H=hermitetx(x0,y0,xi,x,y)H= interp1(x0,y0,xi,'pchip');Hn= interp1(x0,y0,x,'pchip');plot(x0,y0,'o',x,Hn,'-',xi,H,'*',x,y,'-.')legend('节点(xi,yi)', '分段埃尔米特插值函数','插值点(x,H)','被插值函数y')我们也可以直接在在MATLAB 工作窗口编程序，例如，>> x0 =-6:6; y0 =sin(x0); xi = -6:.25:6;yi = interp1(x0,y0,xi,'pchip');x=-6:0.001:6; y=sin(x); plot(x0,y0,'o',xi,yi,x,y,':'), legend('节点(xi,yi)','分段埃尔米特插值函数','被插值函数y=sinx') title(' y=sinx 及其分段埃尔米特插值函数和节点的图形')>> x0 =-6:6; y0 =cos(x0);xi = -6:.25:6;yi = interp1(x0,y0,xi,'pchip');x=-6:0.001:6; y=cos(x); plot(x0,y0,'o',xi,yi,x,y,':'), legend('节点(xi,yi)','分段埃尔米特插值函数','被插值函数y=cosx') title(' y=cosx 及其分段埃尔米特插值函数和节点的图形')例6.6.6 设函数211)(x x f +=定义在区间]5,5[-上，节点(X (i ),f (X (i )))的横坐标向量X 的元素是首项a =-5，末项b =5，公差h =1.5的等差数列，构造三次分段埃尔米特插值函数)(3,x H n .把区间]5,5[-分成20等份，构成20个小区间，用MA TLAB 程序计算各小区间中点i x 处)(3,x H n 的值，并作出节点，插值点，)(x f 和)(3,x H n 的图形.解 在MATLAB 工作窗口输入程序>>x0=-5:1.5:5;y0=1./(1+x0.^2); x1=-4.75:0.5:4.75;x=-5:0.001:5; y=1./(1+x.^2); H= hermitetx(x0,y0,x1,x,y) title('函数y=1/(1+x^2)及其分段埃尔米特插值函数，插值，节点(xi,yi) 的图形')运行后屏幕显示各小区间中点i x 处)(3,x H n 的值，出现节点，插值点，)(x f 和)(3,x H n 的图形(略).例6.6.7 设函数x x x f cos 5.0)(-=定义在区间],[ππ-上，取7=n ，按等距节点构造分段埃尔米特插值函数)(3,7x H ，用MA TLAB 程序计算各小区间中点i x 处)(3,7x H 的值，作出节点，插值点，)(x f 和)(3,7x H 的图形.解 记节点的横坐标7,,2,1,0,7/2, =π=+π-=i h ih x i 插值点)(21121+++=i i i x x x，6,,2,1,0 =i .在MA TLAB 工作窗口输入程序>> h=2*pi/7; x0=-pi:h:pi;y0=0.5.*x0-cos(x0); xi=-pi+h/2:h:pi-h/2;b=max(x0); a=min(x0); x=a:0.001:b;y=0.5.*x-cos(x); H= hermitetx(x0,y0,xi,x,y)title('函数y=0.5x-cos(x)及其分段埃尔米特插值函数，插值，节点(xi,yi) 的图形')运行后屏幕显示各小区间中点i x 处)(3,7x H 的值，出现节点，插值点，)(x f 和)(3,7x H 的图形(略).6.6.4 用MATLAB 计算有关分段埃尔米特插值的误差例6.6.8 设函数22511)(x x f +=定义在区间]1,1[-上，取10=n ，按等距节点构造分段埃尔米特插值函数)(3,x H n ，用MA TLAB 程序在]1,1[-上计算)(max )4(11x fx ≤≤-和)(3,x H n 的误差公式和误差限. 解 在MATLAB 工作窗口输入程序>> syms h,x=-1:0.0001:1;yxxxx=150000000./(1+25.*x.^2).^5.*x.^4-4500000./(1+25.*x.^2).^4.*x.^2+15000./(1+25.*x.^2).^3;myxxxx=max(yxxxx), R=(h^4)* abs(myxxxx/384)运行后输出)(x f 的4阶导数在区间]1,1[-上绝对值的最大值myxxxx 和)(3,x H n 在区间]1,1[-上的误差公式myxxxx 为myxxxx = R =15000 625/16*h^4（4）在MATLAB 工作窗口输入程序>> h=0.2; R =625/16*h^4运行后输出误差限为R =0.06250000000000例6.6.9 设函数))432sin 3tan(cos()(2x x x f ++=定义在区间],[ππ-上，取9=n ，按等距节点构造分段埃尔米特插值函数)(3,x H n .(1)用MA TLAB 程序计算各小区间中点i x 处)(3,x H n 的值，作出节点，插值点，)(x f 和)(3,x H n 的图形； (2) 并用MA TLAB 程序计算各小区间中点处)(3,x H n 的值及其相对误差；(3) 用MA TLAB 程序求)(max )4(x f x ππ≤≤-和)(3,x H n 在区间],[ππ-上的误差公式和各插值的误差限.解 （1）记节点的横坐标9,,2,1,0,9/2, =π=+π-=i h ih x i ，插值点)(21121+++=i i i x x x，8,,2,1,0 =i .在MATLAB 工作窗口输入程序>>h=2*pi/9; x0=-pi:h:pi;y0=tan(cos((sqrt(3)+sin(2*x0))./(3+4*x0.^2)));xi=-pi+h/2:h:pi-h/2;fi=tan(cos((3^(1/2)+sin(2*xi))./(3+4*xi.^2)));b=max(x0); a=min(x0); x=a:0.001:b;y=tan(cos((3^(1/2)+sin(2.*x))./(3+4*x.^2)));Hi= hermite tx (x0,y0,xi,x,y);Ri=abs((fi-Hi)./fi); xi,fi,Hi,Ri,i=[xi',fi',Hi',Ri']title('函数y=tan(cos((sqrt(3)+sin(2x))/(3+4x^2)))及其分段埃尔米特插值函数，插值，节点(xi,yi) 的图形')运行后屏幕显示各小区间中点x i 处的函数值f i ，插值H i ，相对误差值R i ，并且作出节点，插值点，)(x f 和)(3,x H n 的图形(略).（2）在MATLAB 工作窗口输入程序>> syms xy=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2)));yxxxx=diff(y,x,4),%simple(yxxxx)运行后屏幕显示函数)(x f 的4阶导数)()4(x f，然后将输出的)()4(x f 编程求)(m a x )4(x f x ππ≤≤-和)(3,x H n 及其在区间],[ππ-上的误差限的MA TLAB 程序如下>>syms h,x=-pi:0.0001:pi;yxxxx=-12.*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2).^2.*sin((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)).^3.*(2.*cos(2.*x)./(3.+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^2.*(-4.*sin(2.*x)./(3.+4.*x.^2)-32.*cos(2.*x)./(3.+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^3.*x.^2.-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2)+16.*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2).^2.*sin((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)).^4.*(2.*cos(2.*x)./(3.+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^4.*tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)))+8.*tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^3.*(1.+tan(co s((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^4.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^4-8.*tan(cos((3.^(1./2)+sin(2.*x ))./(3.+4.*x.^2))).*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)).^2.*(2.*cos(2.*x)./(3.+4.*x.^2)-8*(3.^(1./2)+sin(2.*x))./(3.+4*x.^2).^2.*x).^4+6.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*(-4.*sin(2.*x)./(3+4.*x .^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)+(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4-3.*(1+tan(cos((3.^(1./2)+si n(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x +128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+s in(2.*x))./(3+4.*x.^2).^2).^2-4.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).*(-8.*cos(2.*x)./(3+4.*x.^2)+96.*sin(2.*x)./(3+4.*x.^2).^2.*x +768.*cos(2.*x)./(3+4.*x.^2).^3.*x.^2-48.*cos(2.*x)./(3+4.*x.^2).^2-3072.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^3+384.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x)-(1+tan(cos((3.^(1./2)+sin(2.*x ))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(16.*sin(2.*x)./(3+4.*x.^2)+256.*cos(2.*x)./(3+4.*x.^2).^2.*x-3072.*sin(2.*x)./(3+4.*x.^2).^3.*x.^2+192.*sin(2.*x)./(3+4.*x.^2).^2-24576.*cos(2.*x)./(3+4.*x.^2).^4.*x.^3+3072.*cos(2.*x)./(3+4.*x.^2).^3.*x+98304.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^5.*x.^4-18432.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^2+384.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3)-12.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).^2*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4.*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))-24.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^3.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x .^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)-24.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4.*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))+36.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)+6.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4+6.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*s in((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2).^2+8.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).*(-8.*cos(2.*x)./(3+4.*x.^2)+96.*sin(2.*x)./(3+4.*x.^2).^2.*x+768.*cos(2.*x)./(3+4.*x.^2).^3.*x.^2-48.*c os(2.*x)./(3+4.*x.^2).^2-3072.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^3+384.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x)myxx=max(yxxxx), R=(h.^4).* abs(myxx./384)将其保存为名为myxx.m 的M 文件，然后在MATLAB 工作窗口输入该文件名>> myxx运行后屏幕显示myxx =)(max )4(x f x π≤≤π-和)(3,x H n 在区间],[ππ-上的误差公式)(max384)4(4x f h R x π≤≤π-=如下myxx = R =73.94706841647552 1734520780029061/9007199254740992*h^4最后在MATLAB 工作窗口输入>>h=2*pi/9; R =1734520780029061/9007199254740992*h^4运行后屏幕显示)(3,x H n 在区间],[ππ-上的误差限R =0.04574453029948。
matlab数值分析第三章插值
• 一个多项式通常不用拉格朗日形式表示，它更 常见的写成类似
x 2x 5
3
• 的形式。其中简单的x的次方项称为单项式， 而多项式的这种形式称为使用幂形式的多项式。 • 插值多项式使用幂形式表示为
P( x) c1x c2 x ... cn1x cn
n1Βιβλιοθήκη n 2• 其中的系数，原则上可以通过求解下面的线性代 数方程组得到。
3.2 分段线性插值
• • • • 通过两步操作可以绘制出一个简单的图形： 第一步用圆圈在坐标系中标出个数据点plot(x,y,'o'); ， 第二步用直线段依次连接这些数据点plot(x,y'-'); 。 下面的语句执行这样的操作，生成图3-3.
• x = 1:6; • y = [16 18 21 17 15 12]; • plot(x,y,'o',x,y,'-');
3.4 保形分段三次插值
• pchip实际是“分段三次埃米特插值多项式”
（piecewise cubic Hermite interpolating polynominal）的
英文首字母缩写。有意思的是，根据这个名字并不能 确定它到底是哪一种分段三次埃米特插值多项式，因 为样条插值函数实际也是分段三次埃米特插值多项式， 只是对斜率的限制条件不同而已。 • 在这里，我们说的pchip实际上是一个最近才引入 MATLAB、保形的（shape-preserving）且看上去不 错的特定插值函数。它基于一个由Fritsch和Carlson 编写的旧的Fortran程序，在Kahaner、Moler和 Nash的书【33】中可以找到相关的介绍。
V=vander(x) 生成 V = 0 0 1 1 8 4 27 9 然后，输入命令 c=V\y' 计算出插值系数 c = 1.0000 0.0000 -2.0000 -5.0000
实习：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插值c++程序例题
两点三次hermite插值c++程序例题两点三次Hermite插值是一种数值分析方法，用于在给定的数据点之间估计函数值。
这种方法基于多项式插值，并使用导数信息来提高插值的准确性。
以下是一个C++程序，实现了两点三次Hermite插值：cpp#include <iostream>#include <vector>// 定义一个结构体，用于存储数据点和它们的导数struct Point {double x;double y;double dy;};// 计算两点之间的差值double difference(double a, double b) {return a - b;}// 计算两点之间的差值的平方double squareDifference(double a, double b) {return difference(a, b) * difference(a, b);}// 计算两点之间的差值的立方double cubeDifference(double a, double b) {return squareDifference(a, b) * difference(a, b);}// 两点三次Hermite插值函数double hermiteInterpolation(const std::vector<Point>& points, double x) {double result = 0.0;for (size_t i = 0; i < points.size(); ++i) {double term = points[i].y;double prod = 1.0;for (size_t j = 0; j < points.size(); ++j) {if (i != j) {double weight = cubeDifference(x, points[j].x) / (squareDifference(points[i].x, points[j].x) *squareDifference(x, points[i].x));prod *= weight;}}result += term * prod;}return result;}int main() {std::vector<Point> points = {{1.0, 2.0, 3.0},{2.0, 3.0, 4.0},{3.0, 5.0, 6.0},{4.0, 7.0, 8.0}};double x = 2.5;double y = hermiteInterpolation(points, x);std::cout << "The interpolated value at x = " << x << " is y = " << y << std::endl;return 0;}这个程序首先定义了一个结构体`Point`，用于存储数据点及其导数。
“Hermie”插值
一、科学计算的算法及其举例应用和利用MATLAB 自带函数实现hermite 插值算法说明：已知n 个插值节点x 1,x 2,…,x n 及其对应的函数值y 1,y 2,...,y n 和一阶导数值y 1',y 2',...,y n '。
则计算插值区域内任意x 的函数值y 的Hermite 插值公式：])2)([()('1i i ni i i i i y y y a x x h x y +--=∑= 其中 ∑∏≠=≠=-=--=n ij i j i i ni j j j i j i x x a x x x x h 1211;)( 源代码： hermite.mfunction 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~=ih=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; a=1/(x0(i)-x0(j))+a; end endyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); end y(k)=yy; end流程图：开始y=hermite(x0,y0,y1,x)n=length(x0) m=length (x) k=1:myy=0.0i=1:nh=1.0 a=0.0j=1:ny(k)=yyyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i))h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2 a=1/(x0(i)-x0(j))+aj~=i结束FalseFalseFalseTrueTrueTrue举例应用：对给定数据表，试构造Hermite多项式，并给出sin0.34的近似值。
x 0.30 0.32 0.35 sinx 0.29552 0.31457 0.34290 cosx 0.95534 0.94924 0.93937流程图：源代码：在MATLAB命令窗中输入>> x0=[0.3 0.32 0.35];>> y0=[0.29552 0.31457 0.34290]; >> y1=[0.95534 0.94924 0.93937]; >> x=[0.3: 0.005: 0.35];>> y=hermite(x0,y0,y1,x);>> plot (x,y)>> y=hermite(x0,y0,y1,0.34)y =输出图形和y值结束plot(x,y) 开始x0=[0.3 0.32 0.35];y0=[0.29552 0.314570.34290]; y1=[0.95534 0.94924 0.93937];x=[0.3: 0.005: 0.35]y=hermite(x0,y0,y1,x)y=hermite(x0,y0,y1,0.34)0.3335>> sin(0.34)ans =0.3335>> y2=sin(x); >> hold on>> plot (x,y2,'--r') 运行结果：二、科学计算和工程实际问题和举例1、求导弹轨迹及加速度流程图：源代码： %构建方程 work1f.mfunction rprime=work1f(t,r) global vt vm rprime=[0;0];% 给出t0之前rprime 初值rprime(1)=-vt-vm*r(1)/sqrt(r(1)^2+r(2)^2); rprime(2)=-vm*r(2)/sqrt(r(1)^2+r(2)^2); %绘制曲线 work1.m global vt vmvt=input('vt=');vm=input('vm='); %输入共用的参数 r0=input('[x0;y0]=');%输入数值积分需要的参数tspan=input('tspan=[t0,tfinal]='); [t,r] = ode45('work1f',tspan,r0); %进行数值积分 plot(r(:,1),r(:,2));%绘图%求加速度：M 点位置的导数是相对速度，二次导数则为绝对加速度dt=diff(t); Ldt=length(dt); %为了求导数先求t 的增量 x=r(:,1);y=r(:,2);% 把r 写成x,y 两个分量形式 vx=diff(r(:,1))./dt;vy=diff(r(:,2))./dt;开始构建方程 数值积分输入全局参数输出图形求2次导数求加速度输出加速度数值解结束wx=diff(vx)./dt(1:Ldt-1);wy=diff(vy)./dt(1:Ldt-1); %二次导数[t(2:Ldt),x(2:Ldt),y(2:Ldt),wx,wy] %显示数据在命令行输入>>work1>>vt=500;vm=1000；>>[x0;y0]=[3000;4000]；>>tspan=[t0,tfinal]=[0,4.5]>>hold on>>work1>>vt=500;vm=800；>>[x0;y0]=[3000;4000]；>>tspan=[t0,tfinal]=[0,6]运行结果：分别获得两组解，获得轨迹图如下两组绝对加速度的数值解按列以对应的时间、相对位置的x、y投影距离和x、y方向绝对加速度输出Vm=800ans =1.0e+003 *0.0002 2.8536 3.9036 0.0539 -0.03940.0003 2.7084 3.8062 0.0569 -0.0405 0.0005 2.5645 3.7080 0.0600 -0.0415 0.0006 2.4219 3.6088 0.0635 -0.0426 0.0008 2.2808 3.5087 0.0672 -0.0437 0.0009 2.1411 3.4076 0.0713 -0.0448 0.0011 2.0031 3.3055 0.0757 -0.0459 0.0012 1.8668 3.2023 0.0806 -0.0469 0.0014 1.7323 3.0981 0.0859 -0.0480 0.0015 1.5997 2.9929 0.0917 -0.0490 0.0017 1.4692 2.8865 0.0982 -0.0499 0.0018 1.3409 2.7790 0.1052 -0.0507 0.0020 1.2149 2.6703 0.1130 -0.0514 0.0021 1.0915 2.5605 0.1217 -0.0518 0.0023 0.9709 2.4495 0.1313 -0.0520 0.0024 0.8532 2.3374 0.1419 -0.0517 0.0026 0.7387 2.2241 0.1538 -0.0510 0.0027 0.6276 2.1096 0.1671 -0.0496 0.0029 0.5203 1.9941 0.1819 -0.0473 0.0030 0.4172 1.8774 0.1983 -0.0439 0.0032 0.3184 1.7598 0.2168 -0.0391 0.0033 0.2246 1.6413 0.2373 -0.0322 0.0035 0.1360 1.5221 0.2599 -0.0228 0.0036 0.0534 1.4023 0.2843 -0.0106 0.0038 -0.0229 1.2824 0.3122 0.0058 0.0039 -0.0921 1.1625 0.3415 0.0277 0.0041 -0.1537 1.0433 0.3711 0.0555 0.0042 -0.2069 0.9253 0.3996 0.0890 0.0044 -0.2511 0.8094 0.4317 0.1348 0.0045 -0.2856 0.6964 0.4570 0.1894 0.0047 -0.3099 0.5877 0.4692 0.2489 0.0048 -0.3235 0.4847 0.4184 0.2734 0.0049 -0.3269 0.4086 0.4678 0.3750 0.0050 -0.3237 0.3378 0.4514 0.4339 0.0052 -0.3144 0.2730 0.4110 0.4739 0.0053 -0.2993 0.2147 0.2959 0.40420.0053 -0.2874 0.1822 0.3379 0.53260.0054 -0.2738 0.1524 0.3056 0.54840.0055 -0.2585 0.1256 0.2653 0.54620.0056 -0.2419 0.1016 0.1991 0.46920.0056 -0.2285 0.0854 0.2022 0.53930.0057 -0.2146 0.0709 0.1770 0.53450.0057 -0.2001 0.0579 0.1483 0.51250.0058 -0.1852 0.0465 0.1062 0.41920.0058 -0.1741 0.0392 0.1085 0.48110.0059 -0.1629 0.0326 0.0940 0.46890.0059 -0.1516 0.0267 0.0783 0.44410.0059 -0.1401 0.0216 0.0458 0.29130.0060 -0.1357 0.0198 0.0606 0.41610.0060 -0.1313 0.0181 0.0562 0.40800.0060 -0.1268 0.0165 0.0518 0.3991Vm=1000ans =1.0e+003 *0.0001 2.8767 3.9097 0.0668 -0.04920.0002 2.7542 3.8188 0.0699 -0.05040.0003 2.6326 3.7272 0.0731 -0.05160.0004 2.5119 3.6350 0.0766 -0.05290.0006 2.3922 3.5421 0.0804 -0.05430.0007 2.2735 3.4485 0.0844 -0.05560.0008 2.1558 3.3542 0.0888 -0.05700.0009 2.0393 3.2592 0.0935 -0.05850.0010 1.9240 3.1635 0.0986 -0.06000.0011 1.8099 3.0670 0.1042 -0.06150.0012 1.6972 2.9697 0.1102 -0.06300.0014 1.5858 2.8716 0.1168 -0.06450.0015 1.4759 2.7727 0.1241 -0.06600.0016 1.3676 2.6730 0.1320 -0.06750.0017 1.2610 2.5724 0.1408 -0.0690 0.0018 1.1561 2.4710 0.1504 -0.0703 0.0019 1.0532 2.3686 0.1611 -0.0716 0.0020 0.9523 2.2654 0.1730 -0.0727 0.0021 0.8535 2.1612 0.1863 -0.0735 0.0022 0.7572 2.0561 0.2011 -0.0739 0.0024 0.6634 1.9501 0.2177 -0.0740 0.0025 0.5723 1.8431 0.2366 -0.0733 0.0026 0.4842 1.7352 0.2578 -0.0717 0.0027 0.3994 1.6264 0.2817 -0.0690 0.0028 0.3182 1.5167 0.3093 -0.0648 0.0029 0.2408 1.4062 0.3411 -0.0581 0.0030 0.1678 1.2950 0.3772 -0.0483 0.0032 0.0996 1.1831 0.4175 -0.0349 0.0033 0.0366 1.0708 0.4670 -0.0157 0.0034 -0.0204 0.9583 0.5238 0.0125 0.0035 -0.0708 0.8460 0.5857 0.0508 0.0036 -0.1138 0.7343 0.6459 0.0949 0.0037 -0.1487 0.6238 0.7458 0.1799 0.0038 -0.1740 0.5156 0.8446 0.2962 0.0039 -0.1887 0.4112 0.9078 0.4245 0.0041 -0.1919 0.3121 0.6878 0.3993 0.0041 -0.1895 0.2699 1.0184 0.7163 0.0042 -0.1844 0.2296 1.0405 0.8403 0.0042 -0.1767 0.1914 1.0333 0.9575 0.0043 -0.1664 0.1556 0.8512 0.8941 0.0043 -0.1579 0.1327 1.0076 1.1994 0.0043 -0.1481 0.1112 0.9850 1.3162 0.0044 -0.1372 0.0913 0.9290 1.3989 0.0044 -0.1252 0.0731 0.7705 1.2960 0.0044 -0.1148 0.0597 0.8337 1.6027 0.0044 -0.1038 0.0476 0.7886 1.7331 0.0045 -0.0923 0.0368 0.6839 1.72432、用周期样条求曲线轮廓并作图流程图：源代码： work2.m function ppx=[100 134 164 180 198 195 186 160 136 100 66 35 15 0 5 17 32 63 100];y=[503 525 514.3 451 326.5 188.6 92.2 59.6 62.2 102.7 147.1 191.6 236 280.5 324.9 369.4 413.8 458.3 503]; i=[0:pi/9:pi*2];xx=csape(i,x,'periodic'); yy=csape(i,y,'periodic'); ii=[0:pi/225:pi*2]; xpp=ppval(xx,ii); ypp=ppval(yy,ii);pp=plot(xpp,ypp,'r-',x,y,'bd'); 运行结果：获得封闭曲线如右图取作图用样点作图输出图形构建样条函数结束开始断点数据。
基于matlab的常见插值法及其应用
基于matlab的常见插值法及其应用郭小乐【摘要】本文就数值分析中几种常见的插值法:拉格朗日插值、牛顿插值、Hermite插值及三次样条插值,讨论其不同形式的表达式及误差,结合matlab给出具体实例,对比分析.此外还就三次样条插值的不同计算方法进行归纳、总结.【期刊名称】《赤峰学院学报（自然科学版）》【年(卷),期】2017(033)007【总页数】3页(P5-7)【关键词】拉格朗日插值;牛顿插值;Hermite插值;三次样条插值;matlab【作者】郭小乐【作者单位】宁夏大学数学统计学院,宁夏银川 750021【正文语种】中文【中图分类】O241.3在许多工程问题中,有时我们只能给出某一函数在一些离散点的值,给不出具体的函数表达式,或者函数的表达式过于复杂不利于计算,这时我们就需要构造这个函数的近似函数,数学上我们把这种方法称为插值[1].插值法作为函数逼近、数值微积分及微分方程数值解的基础,在当今社会被越来越多的学者所关注.尤其随着计算机的普及,很多学者将插值与matlab等软件结合,使得插值法得到了更广泛的应用.常用的插值法包括：拉格朗日 (Lagrange)插值、牛顿(Newton)插值、Hermite 插值、三次样条插值,本文就这四种插值法,结合matlab,从其公式的构造、余项出发,对不同的插值法通过数值实例进行对比研究.此外,还就具有良好收敛性的三次样条函数,归纳出不同的计算方法.由于代数多项式具有简单和一些良好的特性,如多项式是无穷光滑的,容易计算它的导数及积分.故本文选择代数多项式作为插值函数.下面依次就这几种插值法进行讨论.2.1 Lagrange插值所谓n次Lagrange插值即给定平面上的n+1个互不相同的插值点(xi,f(xi)),i=0,1,2,…,n,利用插值基函数构造唯一的一条次数不高于n次的插值多项式.2.1.1 n次Lagrange插值多项式形式及误差拉格朗日插值公式的优点是格式整齐和规范,在理论分析中有重要的价值.它的缺点是没有承袭性质,当需要增加插值节点时,需要重新计算所有的插值基函数,计算繁琐.2.2 Newton插值Newton插值是利用差商计算的,下面我们首先给出差商的概念.定义1[1]设x0,x1,…，xk互不相同,f(x)关于x0,x1,…，xk的k阶差商为2.2.1 n次Newton插值多项式的形式及误差牛顿插值具有拉格朗日插值没有的承袭性,当增加插值节点时,只需再计算一项即可得到相应的插值多项式.由插值多项式的唯一性可知二者是同一插值多项式的不同表达形式.2.3 Hermite插值常用Hermite插值描述如下[1]：对于f(x)具有一阶连续导数,以及插值点xi,i=0,1,…,n,xi互不相同,若有至多为2n+1次的多项式函数H2n+1(x)满足：则称H2n+1(x)为f(x)关于节点{xi}ni=0的Hermite插值多项式.2.3.1 Hermite插值多项式的形式及误差:2.4 三次样条插值在实际计算中,人们很少采用高次插值,其主要原因有两个:其一,由于受到所要通过的插值节点的约束,高次插值描绘的代数曲线是摆动的.即我们常说的Runge(龙格)现象;其二,从计算的舍入误差看,对于等距结点的差分形式,结点上函数值的微小变化可能导致高阶差分的很大变动.故在应用中主要采用的是分段低次插值.实践证明,用分段的低次插值多项式逼近被插函数比在全区间上用高次插值多项式效果好.一般来说,分段插值所描绘的曲线是不光滑的.为了能获得一条足够光滑的插值曲线,近年来人们采用了分段多项式光滑插值法,即著名的样条函数插值法.本文我们主要介绍具有良好的收敛性,具有二阶光滑性的三次样条插值.2.4.1 三次样条函数定义2[5]给定[a,b]上n+1个节点a=x0＜x1＜…＜xn-1＜xn=b以及这些点上的函数值f(xi)=yi(i=0,1,…,n).若函数s(x)满足:（1）s(xi)=yi,i=0,1,2,…,n；（2）在每个小区间[xi,xi+1]上是一个次数不超过三次的多项式；（3）s(x)、s'(x)、s"(x)在[a,b]上都连续.则称s(x)为函数f(x)关于节点x0,x1,…,xn的三次样条函数.2.4.2 三次样条插值的几种计算方法2.4.2.1 利用二阶导数为线性函数求解令s"(xj)=Mj,j=0,1,…,n.因为s(x)为分段三次多项式,故s"(x)在[xj,xj+1]上是线性函数.可表示为其中hj=xj+1-xj为了求出s(x)在[xj,xj+1]上的表达式,需要对上式积分两次,并利用s(xj)=yi及s(xj+1)=yj+1可求出积分常数.得三次样条函数表达式:再利用光滑连接条件s'(xj-0)=s'(xj+0)，并对s(x)求导得为了确定唯一的Mj(j=0,1,…,N),补充三个边界条件:（1）假定s'(x0)=f'0，s'(xn)=f'n；（2）假定 s"(x0)=f"(x0)，s"(xn)=f"(xn)直接得端点方程M0=f"0，Mn=f"n；（3）假定M0=MN=0,即为自然边界条件,可得M0=Mn，λnM1+μnMn-1+2Mn=dn；以上三个边界条件将其写成矩阵形式,得关于Mj(j=0,1,…,n)的三对角线性方程组,利用追赶法[2]求出唯一的Mj(j=0, 1,…,n),将其代入s(x)即可.2.4.2.2 利用特殊形式的埃尔米特插值公式求解若插值节点取为xk,xk+1,插值多项式为H3(x),则采用基函数法,令H3(x)=αk(x)yk+αk+1yk+1+βk(x)mk+βk+1(x)mk+1其中αk(x)，αk+1（x），βk(x)，βk+1(x)是关于插值节点 xk，xk+1的三次Hermite插值基函数.上式即为由特殊的埃尔米特插值公式推导的三次样条插值的计算方法之一.例设，x3=4.利用不同的插值法构造f(x)在[，4]上的插值多项式,使其满足P(xi)=f(xi) (i=0,1,2,3)，P'(x0)=f'(x0).下图是利用matlab(2012b)编程软件绘制的上述条件下四种不同算法在[0.25,4]插值函数曲线图(a)及其相对误差曲线图(b),直观明了地展示不同插值算法的效果.从图a中可看出不同插值方法对原函数的逼近效果相差不大.但从图b中的不同插值法的相对误差大小可知,这四种插值方法中,Lagrange插值和Newton插值的效果是一样的,相对另外两种插值方法来说,这两种方法的逼近效果较差,其次是Hermite插值的逼近效果较好,而具有良好收敛性及二阶光滑性的三次样条插值来说,其误差最小,逼近效果最好.本文根据运用的插值条件的不同,对《数值分析》中常用的四种插值法进行归纳总结.通过具体的数值实例,结合matlab编程软件,通过图像直观明了地展示不同方法的逼近效果,还绘制了不同插值方法的相对误差曲线图,从误差的角度分析不同方法的效果.为不同领域的科研工作者提供理论基础.〔1〕张韵华,奚梅成,陈效群.数值计算方法与算法[M].北京：科学出版社,2006. 〔2〕李庆扬,王超能,易大义.数值分析（第五版）[M].北京：清华大学出版社,2008. 〔3〕徐利治,周蕴时,孙玉柏.逼近论[M].北京：国防工业出版社,1985.〔4〕吕晓亚,张莉.插值法在数值分析中的教学实践[J].唐山师范学院学报,2011,33(2):123-125.〔5〕张希娜,李亚红,郭中凯.关于三次样条插值的教学研究[J].长沙大学学报,2012,26(2):131-132.〔6〕张丽娟.三种插值方法的应用与比较[J].赤峰学院学报(自然科学版),2010,26(3):1-3.〔7〕林昌华,杨岩.拉格朗日插值法在工程设计及CAD中的应用[J].重庆理工大学学报(自然科学版),2013,27(12):34-37.。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
给定矢量P 0, P 1, R 0, R 1,称满足下列条件的参数三次多项式曲线P (t ),t ∈[0,1]为Hermite 曲线：
H (x 0)=y 0,H (x 1)=y 1,
H ′(x 0)=m 0,H ′(x 1)=m 1,
即Hermite 曲线两个端点为P 0, P 1,在两端点的切矢量分别R 0, R 1。
记几何矩阵和基矩阵分别为G H , M H , G H , M H 是未知的.取G H =[P 0,P 1,R 0,R 1],则只要M H 就可以了。
一般的曲线经过多项式分解, 得到参数多项式曲线的矩阵表示:
P (t )=G ∙M ∙T
将（1）式代入（2）得到：
G H ∙M H ∙T H |t=0=G H ∙M H ∙(1,0,0,0)T =P 0, G H ∙M H ∙T H |t=1=G H ∙M H ∙(1,1,1,1)T =P 1, G H ∙M H ∙T H |t=0=G H ∙M H ∙(0,1,0,0)T =R 0, G H ∙M H ∙T H |t=0=G H ∙M H ∙(0,1,2,3)T =R 1, 将上面四个式子合并如下形式：
G H ∙M H ∙[1
01001111020103
]=[P 0,P 1,R 0,R 1]=G H 上面方程的解不唯一，不妨取
M H =[1
01001111020103]−1=[1000−320−3−21−2100−11
] 从而得到三次Hermite 曲线的方程：
P(t)=G H∙M H∙T
其中M H∙T确定了一组Hermite基函数G0(t)，G1(t)，H0(t)，H1(t)，即
M H∙T=[1
0−32
0−3−2
1−21
00−11
][
1
t
t2
t3
]=[
1−3t2+2t3
3t2−2t3
t−2t2+t3
−t2+t3
]
附：MATLAB程序
function yy=hermite(x,y,dy,xx)
% 输入
X——左右两个端点的X轴坐标
Y——左右两个端点的Y轴坐标
dy——左右两个端点的切矢
xx——中间插值的点X轴坐标
%输出
yy——中间插值的点Y轴坐标
function yy=hermite(x,y,dy,xx)
k=length(xx);
z=zeros(1,k);
for i=1:k;
s=0;
xaix=xx(i);
a=1-3.*(xaix)^2+2.*(xaix)^3;
b=2.*(xaix)^2-2.*(xaix)^3;
c=xaix-2.*(xaix)^2+(xaix)^3;
d=-2.*(xaix)^2+(xaix)^3;
s=y(1)*a+y(2)*b+dy(1)*c+dy(2)*d;
z(i)=s;
end
yy=z;
function yy=hermite(x,y,dy,xx) % 输入
X——左右两个端点的X轴坐标
Y——左右两个端点的Y轴坐标
dy——左右两个端点的切矢
xx——中间插值的点X轴坐标
%输出
yy——中间插值的点Y轴坐标
m=length(x);
n=length(y);
l=length(dy);
k=length(xx);
if m~=n,error('向量长度不一样'); end;
if n~=l,error('向量长度不一样'); end;
z=zeros(1,k);
for i=1:k;
s=0;
a=xx(i)-x(1);
b=x(1)-x(2);
c=xx(i)-x(2);
a1=(1-2*a/b)*(c/b)^2;
aa=xx(i)-x(2);
a2=(1+2*aa/b)*(a/b)^2;
b1=a*(c/b)^2;
b2=c*(a/b)^2;
s=y(1)*a1+y(2)*a2+dy(1)*b1+dy(2)*b2; z(i)=s;
end
yy=z;。