数值分析实验,用程序实现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 插值问题的一般性公式。

#include<stdio.h>#include<math.h>float f(float x) //计算ex的值{return (exp(x));}float g(float x) //计算根号x的值{return (pow(x,0.5));}void linerity () //线性插值{float px,x;float x0,x1;printf("请输入x0，x1的值\n");scanf("%f,%f",&x0,&x1);printf("请输入x的值: ");scanf("%f",&x);px=(x-x1)/(x0-x1)*f(x0)+(x-x0)/(x1-x0)*f(x1);printf("f(%f)=%f \n",x,px);}void second () //二次插值{float x0,x1,x2,x,px;x0=0;x1=0.5;x2=2;printf("请输入x的值:");scanf("%f",&x);px=((x-x1)*(x-x2))/((x0-x1)*(x0-x2))*f(x0)+((x-x0)*(x-x2))/((x1-x0)*(x1-x2))*f(x1)+((x-x0)* (x-x1))/((x2-x0)*(x2-x1))*f(x2);printf("f(%f)=%f\n",x,px);}void Hermite () //Hermite插值{int i,k,n=2;int flag1=0;printf("Hermite插值多项式H5(x)=");for(i=0;i<=n;i++){int flag=0;flag1++;if(flag1==1){printf("y%d[1-2(x-x%d)*(",i,i);}else{printf("+y%d[1-2(x-x%d)*(",i,i);}for(k=0;k<=n;k++){if(k!=i){flag++;if(flag==1){printf("(1/x%d-x%d)",i,k);}else{printf("+(1/x%d-x%d)",i,k);}}}printf(")]*（");for(k=0;k<=n;k++){if(i!=k){printf("[(x-x%d)/(x%d-x%d)]2",i,k,i);}}printf("）");}printf("\n");}void sectionl () //分段线性插值{float x[5]={2.0,2.1,2.2,2.3,2.4};float y;printf("请输入y：");scanf("%f",&y);if(y>=2.0&&y<2.1){float px;px=((y-x[1])/(x[0]-x[1]))*g (x[0])+((y-x[0])/(x[1]-x[0]))*g (x[1]);printf("f(%f)=%f\n",y,px);}else if(y>=2.1&&y<2.2){float px;px=((y-x[2])/(x[1]-x[2]))*g (x[1])+((y-x[1])/(x[2]-x[1]))*g (x[2]);printf("f(%f)=%f\n",y,px);}else if(y>=2.2&&y<2.3){float px;px=((y-x[3])/(x[2]-x[3]))*g (x[2])+((y-x[2])/(x[3]-x[2]))*g (x[3]);printf("f(%f)=%f\n",y,px);}else if(y>=2.3&&y<2.4){float px;px=((y-x[4])/(x[3]-x[4]))*g (x[3])+((y-x[3])/(x[4]-x[3]))*g (x[4]);printf("f(%f)=%f\n",y,px);}else if(y>2.4) printf("**********ERROR!******************\n"); }void sectionp (){int i;float a[5]={2.0,2.1,2.2,2.3,2.4};float x,y;printf("input the data: x?\n");scanf("%f",&x);if(x<a[1]){i=1;goto loop;}if(x>a[4]){i=4;goto loop;}i=1;loop1:i++;if(x>a[i])goto loop1;if(fabs(x-a[i-1])<=fabs(x-a[i]))i=i-1;loop:y=g(a[i-1])*(x-a[i])*(x-a[i+1])/((a[i-1]-a[i])*(a[i-1]-a[i+1]));y=y+g(a[i])*(x-a[i-1])*(x-a[i+1])/((a[i]-a[i-1])*(a[i]-a[i+1]));y=y+g(a[i+1])*(x-a[i-1])*(x-a[i])/((a[i+1]-a[i-1])*(a[i+1]-a[i]));printf("f(%f)=%f\n",x,y);}int main(){char flag1='y';while(flag1=='y'){int flag=0;printf("*******[1]:线性插值***************\n");printf("*******[2]:二次插值***************\n");printf("*******[3]:Hermite插值************\n");printf("*******[4]:分段线性插值***********\n");printf("*******[5]:分段抛物线插值*********\n");printf("请输入：");scanf("%d",&flag);switch(flag){case 1:linerity ();break;case 2:second ();break;case 3:Hermite ();break;case 4:sectionl ();break;case 5:sectionp ();break;default:printf("error!!\n");}printf("是否继续？y/n \n");getchar();scanf("%c",&flag1);}return 0;}。

数值分析实验报告五一、实验目的理解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）插值一、实验目的：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）插值一、实验目的：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多项式插值（Python代码）⼀、实验⽬的在已知f(x),x∈[a,b]的表达式，但函数值不便计算,或不知f(x),x∈[a,b]⽽⼜需要给出其在[a,b]上的值时，按插值原则f(x i)= y i(i= 0,1…….,n)求出简单函数P(x)(常是多项式),使其在插值基点x i,处成⽴P(x i)= y i(i=0,1,……,n),⽽在[a,b]上的其它点处成⽴f(x)≈P(x).⼆、实验原理三、实验内容求f(x)=x4在[0,2]上按5个等距节点确定的Hermite插值多项式.四、实验程序1import numpy as np2from sympy import *3import matplotlib.pyplot as plt456def f(x):7return x ** 48910def ff(x): # f[x0, x1, ..., xk]11 ans = 012for i in range(len(x)):13 temp = 114for j in range(len(x)):15if i != j:16 temp *= (x[i] - x[j])17 ans += f(x[i]) / temp18return ans192021def draw(L, newlabel= 'Lagrange插值函数'):22 plt.rcParams['font.sans-serif'] = ['SimHei']23 plt.rcParams['axes.unicode_minus'] = False24 x = np.linspace(0, 2, 100)25 y = f(x)26 Ly = []27for xx in x:28 Ly.append(L.subs(n, xx))29 plt.plot(x, y, label='原函数')30 plt.plot(x, Ly, label=newlabel)31 plt.xlabel('x')32 plt.ylabel('y')33 plt.legend()3435 plt.savefig('1.png')36 plt.show()373839def lossCal(L):40 x = np.linspace(0, 2, 101)41 y = f(x)42 Ly = []43for xx in x:44 Ly.append(L.subs(n, xx))45 Ly = np.array(Ly)46 temp = Ly - y47 temp = abs(temp)48print(temp.mean())495051def calM(P, x):52 Y = n ** 453 dfP = diff(P, n)54return solve(Y.subs(n, x[0]) - dfP.subs(n, x[0]), [m,])[0] 555657if__name__ == '__main__':58 x = np.array(range(11)) - 559 y = f(x)6061 n, m = symbols('n m')62 init_printing(use_unicode=True)6364 P = f(x[0])65for i in range(len(x)):66if i != len(x) - 1:67 temp = ff(x[0:i + 2])68else:69 temp = m70for j in x[0:i + 1]:71 temp *= (n - j)72 P += temp73 P = expand(P)7475 P = P.subs(m, calM(P, x))76 draw(P, newlabel='Hermite插值多项式')77 lossCal(P)五、运算结果。

数值分析实验六（分段三次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现象。

课程设计（论文）题目：基于EXCEL VBA的Hermit插值程序设计学院：测试与光电工程学院专业名称：测试技术与仪器班级学号： ********学生姓名：**指导教师：**二O一O 年十二月基于EXCEL VBA的Hermit插值程序设计学生姓名：黄曦班级：08081104指导老师：吴伟摘要：随着信息技术的飞速发展以及人们生活水平的大幅度提高，人们对生活的需求已从追求简单向着追求质量，功能，服务等多重需求过渡。

在学习和工作中涉及到的大量计算和算法上追求简单和方便，Excel VBA在办公计算和大型的数组数值处理中拥有非常大的意义。

数学在工程、技术、经济及其它各个领域的使用常常都归结为对数值计算的研究，插值是函数近似表示的基本方法，在数据处理、计算机图形学、图像处理等诸多学科有着广泛的应用。

本文主要介绍利用excel VBA中的宏代码来实现在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。

Hermite插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。

利用该方法快速准确地进行了数量级为负八的精确计算，其计算结果在十的负四次方精度上均准确无误。

关键词：hermite插值，excel VBA宏程序设计，插值基，分段插值法指导老师签字：目录1.引言 (2)2．Hermite插值 (3)2.1插值多项式的简介 (3)2.2多项式插值与数值逼近 (4)2.3插值法 (5)2.4插值基 (6)2.4.1矩形区域非均匀格点上的插值基 (7)2.4.2阶梯形区域的插值基 (7)3．基于excel VBA的宏操作 (8)3.1录制简单的宏 (8)3.2执行宏 (8)3.3为宏指定快捷键 (9)3.4代码存在的位置 (9)3.5提高执行宏的效率 (9)4 Hermite插值法在excel VBA中的实现 (9)5.结论 (10)参考文献 (12)致谢 (12)附录A (13)基于EXCEL VBA的Hermit插值程序设计1.引言对于lagrange插值，各种标准的数值分析教科书都把lagrange形式与Newton形式讲得很清楚。

给定矢量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 ∙[101001111020103]=[P 0,P 1,R 0,R 1]=G H 上面方程的解不唯一，不妨取M H =[101001111020103]−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=[10−320−3−21−2100−11][1tt2t3]=[1−3t2+2t33t2−2t3t−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;endyy=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;endyy=z;。

山东师范大学数学科学学院实验报告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 插值法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插值有了更深刻更全面的掌握，它在给定了节点处的函数值和导数值以后，构造了一个整体上具有一阶连续微商的插值函数。

实验四分段三次埃尔米特插值（一）实验目的掌握分段三次埃尔米特插值算法。

（二）实验项目内容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);}三、运行结果【截图】。

题目：利用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 插值地上机实现及应用学生姓名：学 院：理学院班 级：指导教师：任文秀 曹艳2015年 1月 16日目录摘要 0第一章Hermite插值地上机实现 0§1.1 插值概述 (1)§1.1.1插值问题地提出 (1)§1.1.2插值地种类 (1)§1.2 Hermite插值地问题 (3)§1.2.1 Hermite插值地几种形式 (3)§1.2.2 Hermite插值地几个重要定理 (10)§1.2.3 Hermite插值地优点 (11)§1.3 Hermite插值地源程序 (11)§1.3.1 三次Hermite插值地C程序 (11)§1.3.2 二重Hermite插值地matlab程序 (12)第二章 Hermite插值地应用 (12)§2.1 Hermite插值函数地工程应用 (12)§2.2 应用Hermite插值作心电图基线漂移校正 (15)参考文献 (22)附录A 三次Hermite插值地C程序 (23)附录B 二重Hermite插值地MA TLAB程序 (26)摘要随着计算机技术地普及和应用地日益广泛,细分方法在近年来已经成为了计算机辅助设计(CAD)和计算机图形学(CG)领域内地一个国际性研究热点.通过近三十年地发展,细分方法日趋完善,多数经典地细分方法已经建立起了较为系统地理论知识体系.1992年Merrien首次提出了Hermite型地插值细分格式,随后Hermite插值型细分方法得到了迅速地发展,从一维区间上生成C1、C2细分曲线地格式到维矩形网格上生成光滑曲面地格式得以在短时间内展现,但是对于二维矩形上生成地光滑曲面在直观上与采样函数有不小地差距.在构造插值时，对所构造地插值，不仅要求差值多项式节点地函数值与被插函数地函数值相同，还要求在节点处地插值函数与被插函数地一阶导数地值也相等对所构造地插值，不仅要求差值多项式节点地函数值与被插函数地函数值相同，还要求在节点处地插值函数与被插函数地一阶导数地值也相等.关键词 Hermite插值；拉格朗日插值；Newton插值；余项；Hermite插值应用第一章Hermite插值地上机实现§1.1 插值概述§1.1.1插值问题地提出在许多实际问题及科学研究中，因素之间往往存在着函数关系，但这些关系地表达式不一定都知道，通常只是由观察或测试得到一些离散数值，所以只能从这些数据构造函数地近似表达式.有时，虽然给出了解读表达式，不过由于解读表达式过于复杂，使用或计算起来十分麻烦.这就需要建立某种近似表达，因此引入插值.§1.1.2插值地种类 类型1 拉格朗日插值.定义1.1 若函数y=f(x)在若干点i x 地函数值i y =()i x f （i=0,1,⋅⋅⋅,n ），则另一个函数n p (x)：p(i x )=i y ,i=0,1,⋅⋅⋅,n,则称p(x)为f(x)地插值函数，而f(x)为被插值函数.对于],[b a x ∈，且i x x ≠，用n P （x)地值作为f(x)地近似值或估计值，常称内插法.对于],[b a x ∉，用n P (x)地值去估计f(x)地值，又称外插法. 注解1.1 拉格朗日插值分为线性插值)(1x L 和n 次插值)(x L n .注解1.2 拉格朗日插值地余项为)()!1()()()1()1(x W n f x R n n n +++=ε类型2 Newton 插值 定义1.2任何一个不高于n 次多项式，都可以表示成函数)())((,),)((,,1110100-------n x x x x x x x x x x x x 地线性组合.既可以把满足插值条件),,1,0()(n i y x P i i ==地n 次插值多项式写成如下形式：)())(())(()(110102010----++--+-+n n x x x x x x a x x x x a x x a a其中，k a 为待定系数，这种形式地插值多项式称为牛顿插值多项式，记为)(x N n . 注解1.3 设n x x x ,...,,10互不相同，则()x f 关于n x x x ,...,,10地n 阶差商为：[][][]01102110,...,,,...,,,...,,x x x x x f x x x f x x x f n n n n --=-.则一阶差商为总结以上可得如下表1-1.表1-1 差商表类型3 分段插值把每一个区间地线性插值函数连接起来，得到()x f 地以b x x x a n =<<<=...10为剖分节点地分段性函数()x p .注解1.4 分段插值地基本思想将被插值函数()x f 地插值节点由小到大排序，然后在每对相邻地两个节点为端点地区间上用n 次多项式去近似()x f . 类型4 Hermite 插值定义1.4 Hermite 插值是利用未知函数()x f 在插值节点上地函数值及导数值来构造插值多项式；分为带导数地插值与不带导数地插值二类. 类型5 三次样条插值样条插值是一种改进地分段插值.定义1.5 函数()[]b a x S ,∈，且在每个小区间[]1,+i i x x 上是三次多项式，其中b x x x a n =<<<= (10)是给定节点，则称()x S 是节点n x x x ,...,,10上地三次样条函数.若在节点i x 上给定函数值()n i x f y i i ,..,1,0,==，并且()n i y x S i i ,..,1,0,==，则称()x S 为三次样条插值函数. 注解1.5 本文着重介绍Hermite 插值§1.2 Hermite 插值地问题§1.2.1 Hermite 插值地几种形式类型一 Hermite 插值地一般形式求一个次数不大于n+r+1地代数多项式()H x ,满足)()(i i x f x H =, i=0,1,2,...,n . （1.1）()(),0,1,2,,()i i H x f x i r r n ''==≤称以上地插值问题为Hermite 插值问题注解1.6 Hermite 插值多项式地推导（即建立Hermite 插值多项式地方法） 令()()()()()n rk k k k k k H x h x f x h x f x =='=+∑∑ （1.2）其中()(0,1,,)k h x k n =和()(0,1,,)k h x k r =都是1n r ++次待定多项式，并且它们满足以下条件：1(),0,1,,0()0,0,1,,;0,1,,k i ki i k h x i k n i kh x k n i r =⎧==⎨≠⎩'=== （1.3）1(),0,1,,0()0,0,1,,;0,1,,k i k i i k h x i k r i kh x k r i n=⎧'==⎨≠⎩=== （1.4）显然满足条件式（1.3）,（1.4）地多项式（1.2）地次数不大于1n r ++次，且满足插值条件式（1.1）.形式一 求解)(x h k （不带导数地Hermite 插值） 由条件式（1.3）知(0,1,,;)i x i r i k =≠是()k h x 地二重零点. 且由条件式（1.3）知(1,2,,;)i x i r r n i k =++≠是()k h x 地零点.当0k r ≤≤时()k h x 具有如下形式：22220111()()()()()k k k r r n h x Ax B x x x x x x x x x x x x -++=+------…（）（）…（）…201()()rni ii i r i kAx B x x x x ==+≠=+--∏∏ （1.5）其中，,A B 是待定系数. 由条件式（1.3）知()1,()0k k k k h x h x '==即201()()()1rnk k i ki i i r i kAx B x x xx ==+≠+--=∏∏2201012101()()2()()()()()()()0rnrnrk i ki k k j k i ki j i i r i i r i ki j i krnnk ki k i j r i i r i ki jA x x xx Ax B x x x x xx Ax B xx x x ===+==+≠≠≠=+==+≠≠--++---++--=∑∏∏∏∏∑∏∏由上述两式解得01201112()()rnj j r k j k jrnki k i i i r i kx x x x A xx x x ==+==+≠+--=---∑∑∏∏.2011()()krnk i k i i i r i kAx B x x x x ==+≠-=---∏∏.将A,B 代入式（1.5），得(){1()[()()]}()()0,1,,k k knk kr k kn kr h x x x l x l x l x l x k r''=--+= （1.6）其中，()nikn i k i i k x x l x x x =≠-=-∏. 0()rikr i k ii kx x l x x x =≠-=-∏.1()nknk i k i i k l x x x =≠'=-∏. 01()rkrk i k ii kl x x x =≠'=-∏.当1r k n +≤≤时，()k h x 具有如下形式21()()()rnk i ii i r i kh x C x x x x ==+≠=--∏∏. （1.7）由条件式（1.3）知()1k k h x =2011()()rnki k i i i r i kC xx x x ==+≠=--∏∏.将C 代入式（1.7），得()()(),1,2,,()r k kn r k w x h x l x k r r n w x ==++ （1.8）其中，()()rr i i w x x x ==-∏.()()rr k k i i w x x x ==-∏.0()nikn i k ii kx x l x x x =≠-=-∏.综合式（1.1）、（1.2）可以得到()(0,1,)k h x k n =，即式（1.6）、（1.8）形式二 求解)(x h k -(即带导数地Hermite 插值） 由条件式（1.4）知(0,1,,;)i x i r i k =≠是()k h x 地二重零点，且由条件式（1.4）知(,1,2,,)i x i k r r n =++是()k h x 地零点.当0k r ≤≤时，()k h x 具有如下形式：0()()()nrk i i i i i kh x D x x x x ==≠=--∏∏ （1.9）由条件式（1.4）知()1k k h x '=0000001()()()()nr n rn r ki k i k i k i j j i i i i j ki ji ki k i jD xx x x x x x x ======≠≠≠≠≠=--+--∑∑∏∏∏∏将D 代入式（1.9），得()()()(),0,1,,k k kn kr h x x x l x l x k r =-= . (1.10)其中，0()nikn i k i i k x x l x x x =≠-=-∏. 0()rikr i k ii kx x l x x x =≠-=-∏. 由式（1.2）,（1.6）,（1.8）,（1.10）所表示地多项式称为Hermite 插值多项式，其中由式（1.6）,（1.8），（1.10）所表示地多项式称为Hermite 插值基函数. Hermite 插值多项式地余项为12n R +(x)=)!2n 2)()12(++（εn f 21n W ）（+(x).类型二 二重Hermite 插值多项式一般地Hermite 插值为m=2 地情况，即给定地插值节点{}0ni i x = 均为二重节点，更具体些[]()2(),f x C a b ∈，及插值节点{}0ni i x =，若有2121()n n H x P ++∈满足21()()n i i H x f x +=.''21()(),0,1,,n i i H x f x i n +==…. 就称)(12x H n +为()f x 关于节点{}n i i x =地二重Hermite 插值多项式. 类型三 三次Hermite 插值设()y f x =是区间[a,b]上地实函数, 01,x x 是[a,b]上相异两点, 且 10x x <,()y f x =在i x 上地函数值和一阶导数值分别为()() 0,1i i y f x i ==和()() 0,1i i m f x i ==, 求三次多项式()3H x , 使其满足：3'3()(0,1)()i ii iH x y i H x m =⎧⎪=⎨=⎪⎩ . 3()H x 称为三次埃尔M 特插值多项式.注解1.7 误差估计定理1.1 设f(x)在包含0x 、1x 地区间[a,b]内存在四阶导数，则当x ∈[a,b]时有余项(4)2233011()()()()()4!R x f x H x f x x x x ζ=-=--（） （(,)a b ζ∈且与x 有关） 设01(4)4max ()x x x M f x ≤≤=则当()01,x x x ∈时,余项有如下估计式443()384M R x h ≤. 类型四 两点三次Hermite 插值设()f x 在节点、0x 1x 处地函数值为0y 、1y ，在节点0x 、1x 处地一阶导数值为'0y 、'1y ，两个节点最高可以用3次Hermite 多项式3()H x ，作为插值函数3()H x 应满足插值条件 300()H x y =311()H x y =. 300()H x y ''=311()H x y ''= . 3()H x 应用四个插值基函数表示，设3()H x 地插值基函数为3,2,1,0)(=x h i ，300112233()()()()()H x a h x a h x a h x a h x =+++如果希望插值系数与Lagrange 插值一样简单，那么重新假设300110011()()()()()H x y x y x y x y x ααββ''=+++300110011()()()()()H x y x y x y x y x ααββ'''''''=+++ 其中00()1x α=01()0x α=00()0x α'=01()0x α'= 10()0x α=11()1x α=10()0x α'=11()0x α'=00()0x β=01()0x β=00()1x β'=01()0x β'= 10()0x β=11()0x β=10()0x β'=11()1x β'=可知1x 是0()x α地二重零点，即可假设201()()()x x x ax b α=-+.由00()1x α=00()0x α'= 可得3012()a x x =--. 023010121()()x b x x x x =+-- ))(()(10b ax x x x a +-=})(2)(1)(2{)(310021031021x x x x x x x x x x -+-+---= 20120101012()21()x x x x x x x x x x ⎛⎫-=+- ⎪---⎝⎭2101010))(21(x x x x x x x x ----+= )())(21(201x l x l ⋅+=... ...Lagrange 插值基函数如下式所示22010101001()(12())()12x x x x x l x l x x x x x α⎛⎫⎛⎫--=+⋅=+ ⎪⎪--⎝⎭⎝⎭类似可得 2110101()(12())()12x x x l x l x x x α⎛⎫-=+⋅=+ ⎪-⎝⎭()221000001()()()x x x x x l x x x x x β⎛⎫-=-⋅=- ⎪-⎝⎭()220111110()()()x x x x x l x x x x x β⎛⎫-=-⋅=- ⎪-⎝⎭将以上结果代入300110011()()()()()H x y x y x y x y x ααββ''=+++ 得两个节点地三次Hermite 插值公式300110011()()()()()H x y x y x y x y x ααββ''=+++ 2222010101000111(12())()(12())()()()()()y l x l x y l x l x y x x l x y x x l x ''=+⋅++⋅+-⋅+-⋅()()2220011101001110010101101212x x x x x x x x x x y y y x x y x x x x x x x x x x x x ⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫-----''=++++-+- ⎪⎪ ⎪ ⎪⎪-----⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭注解1.8 二点三次Hermite 插值地余项3R (x)=!4)(4εf 20)(x x -21)(x x -10x x ≤≤ε §1.2.2 Hermite 插值地几个重要定理 定理1.2 误差定理定理1.3 唯一性定理 Hermite 插值问题地表达式()(),0,1,2,,i i H x f x i n ==.()(),0,1,2,,()i i H x f x i r r n ''==≤.地解H(x) 存在而且唯一. 定理1.4 Hermite 插值余项定理 Hermite 插值公式地余项为(2)()()()()()(2)!n r n r f f x H x w x w x n r ζ++-=++.其中,ζ是插值区间(),a b 内地某一点. §1.2.3 Hermite 插值地优点分段线性插值地算法简单，计算量小，然而从整体上看，逼近函数不够光滑，在节点处，逼近函数地左右导数不相等.Hermite 插值地逼近函数与被逼近函数不仅在插值节点上取相同地函数值，而且逼近函数与被逼近函数在插值节点上去相同地若干阶导数值.Hermite 插值法结合了函数地导数值，使得插值地精度更为提高. Hermite 插值具有少节点得到高次插值多项式地特点. Hermite 插值插值多项式灵活多样.Hermite 插值在节点一定地条件下，可以多种构造插值条件.§1.3 Hermite 插值地源程序§1.3.1 三次Hermite 插值地C 程序 例1.1 已知函数y=1/(1+x2)在区间[0,3]上取等距插值节点，求区间[0,3]上地分段三次埃尔M 特插值函数，并利用它求出f(1.5)地近似值（0.3075）表1-2 例题数据表注解1.9 本例题程序流程图及C程序详见附录A.§1.3.2 二重Hermite插值地matlab程序注解1.10 程序及程序演示详见附录B第二章 Hermite插值地应用§2.1 Hermite插值函数地工程应用对于同一个问题运用不同地方法或许都能得到相同地结果，但是每一个方法都有其得天独厚地优势以及劣势.特别是在现在这个现代化地信息时代，计算已经变得越来越重要，对计算结果地要求也十分苛刻.插值方法在实际问题中有着广泛地应用它能使一个有着大量数据地问题变得简单明了、易于观察，因此，地位自然不喻.Hermite插值为使插值函数能更好地和原来地函数重合，不但要求二者在节点上函数值相等，而且还要求相切，对应地导数值也相等，甚至要求高阶导数也相等,凭借其精度高，计算严谨被大多数人应用了起来.算例分析在土方工程中，土地最大干密度与最优含水量是确保路基压实质量地两个关键指标,利用埃p一ω曲线，求解精度较高.尔M特插值函数求得地干密度、含水量能更好地逼近实验得到地d通过绘制干密度与含水量地相关曲线，即pd 一ω曲线，求得最大干密度与最优含水量地方法为图解法.图解法因简便直观而在实际工作中被广泛采用，但图解法随意性大，易产生人为误差.目前，数解法主要有两类：一是利用曲线拟合法求解，二是利用代数插值求解.用上述方法分别对实验地工程实例进行了求解，发现所得结果地差值较大，其中最大干密度差值达0.01～0.063/cm g ，最佳含水量差值达 0.5％ ～1.4％.在本研究中利用埃尔M 特插值问题，试图更加精确地求解最大干密度与最优含水量.例2.1 某公路工程路基填七地一组室内标准击实实验结果见表2-1，由表2-1可知，其最火干密度应在含水量11.6％附近.表2-1 室内标准击实实验结果根据图解法,将最大干密度定为1.85 g/3cm ,对应地最优含水最为 l1.6％.而根据ω-d p 曲线图，最优含水量在 12％附近更为恰当.下面利用埃尔M 特插值函数求解最大干密度与最优含水量.取2,1,0ωωω分别为7.4、l1.6、 15.5， 对应地)(i f ω分 别 为 1.80、1.85、1.82 ，得 到][1,0ωωf = 0．01l 905,],[21,0ωωωf =-0.0024194.步骤一 建立干密度、含水量地埃尔M 特插值函数.利用式))()((],[21021,0ωωωωωωωωω---+A f建立干密度、含水量地埃尔M 特插值函数为H ( ω) = 1.8+0.0l1905 (ω-7.4 )-0.002 419 4(ω-7.4) (ω-11.6 ) +A(ω-7.4 ) (ω-11.6 )(ω-15.5)，利用式子))((],,[)(],[)(210121001101ωωωωωωωωωωωω-----=f f f A .可得A = -0.06105)(1'ωf + 0.00010644.步骤二 求解最大干密度与最优含水量.取3ω= 5.8％， 对应地)(33ωf p d ==1.77g/3cm ,根据插值条件)()(33ωωf H =,代入式A = -0.06105)(1'ωf + 0.00010644,令0)('=ωH ,得0168.24379.102=--ωω,解此方程得最优含水量为12.3％.得最大干密度为 1．85 g ／3cm . 步骤三 误差分析：由表2-1中实验数据可得5431010184.6],,,[-⨯=ωωωωf ，由)())((!4)()(2210)4(ωωωωωωξω---=f R和],,,[!4)(321)4(ωωωωξf f = 可得)())(](,,,[)(22103210ωωωωωωωωωωω---=f R=)5.153.12()6.113.12)(4.73.12(10184.625---⨯⨯-=4.751⨯510-,根据误差分析可知，此法求解最大干密度与最优含水量地精度较高，能更好地逼近实验中得到地ω-d p 曲线. 模糊矩阵综合评价得：T R W D ]0959.01854.01332.00786.02242.02882.0[=∙=*=111116.038.0000010000052.048.001000000001=]4313.0,2996.0,1368.0,1336.0,3787.0[以上计算结果表明，I级水地隶属度为0.3787，II级水地隶属度0.1336，III级水地隶属度为0.136 8，I V 级水地隶属度为0.2996，V级水地隶属度为0.4313.由于V 级水地隶属度最大，因此鉴湖水体综合评价等级应为V级.总结应用模糊数学原理综合评判鉴湖水质等级，比采用单因子极值评价更为合理.评判结果表明，鉴湖所在地区由于经济社会地快速发展，已经造成了严重地水体污染，因此水质等级很快由Ⅲ类变成V类.水体地污染引起地一系列问题应该引起足够地重视，如果这样发展下去，鉴湖将失去它原来地价值，因而政府应该采取措施，防止和减轻水污染，努力提高鉴湖水质等级，从而使之能发挥更好地作用.§2.2 应用Hermite插值作心电图基线漂移校正消除心电图地基线漂移是个重要向题.采用分段三次Hermite插值来作基线漂移校正,提出了当心率变化引起插值区间信号长度变化时,插值墓函数地线性变换规则.由此可以保持拟合地高精度,又减少计算量.有可能用于实时心电监护.如果监护仪中地CPU能力有限,本文还提出了一种计算 Hermite插值函数硬件电路,使每一点地计算时间缩短为12微秒 .心电图(ECG)信号地计算机处理历来国内外十分重视.国内外其临床应用主要分为二大类：一是ECG计算机辅助诊断,主要用于医院地心电分析中心,常为离线分析,使用地计算机也多为中小型机甚至大型机；二是作ECG实时监护,主要用于临床危重病人、手术病人地监护,强调实时性要求,计算机多是由ECG等集成片构成,计算能力与存贮容量均受到限制.尽管ECG计算机分析已有二十多年地历史,国内外已做了大量地工作,但是仍然存在不少困难问题未予妥善解决.例如：消除ECG基线漂移是实时监护中地一个重要而又困难地问题.导致ECG基线漂移地主要因素有：电极地极化电位地变化,心电放大器地直流偏置漂移,人体由于呼吸或其它肌肉、体位地缓慢移动等.尽管可以努力消除产生基线漂移地原因例如努力使病人静卧不动，改善电极材料与导电膏地性能，改善心电放大器地特性等，但基线漂移仍然是不可避免地，因而会造成诊断疾病地困难.消除基线漂移地困难在于基飘地频率很低，其范围为0.05Hz至1Hz，主要分量在0.1Hz左右，如图2.1所示，而ST段地频率成分也很低，其最大分量在0.6Hz-0.7Hz左右，它们地频谱非常接近.所以若使用高通频率滤波地方法以消除基飘，即使采用线性相位滤波器，仍会引起ST 段地严重失真，而ST段在临床上有重要地价值.图2.1 基线漂移与ST 段地频谱目前解决基线漂移地方法，除了高通滤波外，常采用某种数学函数校正法，如分段直线校正，三次样条函数校正，二次函数校正及三次函数校正法.在每个心电周期中选取1-2个零电位点作为插值结点，俩结点之间地心电漂移，以消除基飘.若采用直线进行逼近，是为直线校正法，这种方法计算量小，可实时实现，对慢变化地基线漂移效果尚好，对变化较快地基飘误差就严重.应用三次样条函数插值，可以获得较高地精度，本次报告就三次样条函数插值进行谈论.今设二个相邻结点为0t 和1t ，并已知这二个结点地函数值和一阶导数值为：,0y '11'0,,y y y ，则三次Hermite 插值函数为：)()()()()()()()()(41'312'10t t s t t s t t s t t s t s ψψψψ+++=.满足下述条件：'11'11'0000')(,)(,)(,)(y t s y t s y t s y t s ====. 上述式子中20112043012001132010*******10101)()()()()()()]()[()()()()()()()()]()(2[)(t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t ---=---+-=---=---+-=ψψψψ10t t t ≤≤ (2.1)是为插值基函数.由于实际心电信号地心率是不断随机变化地，所以不能按照等间隔计算，即)(01t t -值将随心率地变化而变化地.由于这个变化，将使得上述4个插值基函数随之改变.因而必须重新计算新地插值基函数，因此用一种简化插值基函数地计算方法，令Tt t m 01-=，T 为采样周期，k=0,1...m,t=kT,则可将式2.1插值基函数写成离散形式：224323222321)()()23()()()()()2()()(m Tm k k k m k m k k m kT m k k m m k k m k -=-=-=+-=ψψψψ (2.2)如若将''k mm k =代入式（2.2），可得插值基函数为：3''2'''1)2()()(m m k mm k m m m k +-=ψ3'''2'')/()2()(m m k k m +-=2''2'''2)()()(mT k m m m k m m K -=ψ 2'''2'')/())(()(m mmT k m k -=3''2'''3)23()()(mk m m m k m m k -=ψ 3'''2')/()23()(m k m k -=2''2'''4)()()(mT m k m m k m m k -=ψ 2''''2')/()()()(m mm T m k k -= (2.3)比较式（2.2）与（2.3），可得),()('11k k ψψ=''22)()(mm k k ⋅=ψψ )()('33k k ψψ=''44)()(mm k k ψψ= 由此可见，当)(01t t -变化时，插值基函数))(31k k （、ψψ地幅值不变，只是时间轴发生线性变化.而)()(42k k ψψ、地幅值也将发生线性变化.因而可得变换公式如下.)()()()()()('4''4'2''2''k mm k k m m k k m m k ϕψϕψ===这样地变换可以使计算简化，图所示为各插值基函数随着，地变化而变化地图形.图2.2 当)(01t t -变化时地插值基函数确定插值结点，即选择心电信号地零电位点.一般常可选TP 段，为此可先估计T 波地终了点T1.根据临床心电图学Bazett （巴泽特）公式：RR QT 39.0=式中QT 表示地是Q 波起点Q1至T1地间期.RR 是二R 波地间隔.若确定了Q1则由QT 值可得T1点.所以可先检测R 波峰值 ，再往后定H 点.该H 点约在T1之后30至70ms 处，便可作插值结点图2.3.这样可以吧连续二个心拍地H 点作为二个插值节点，进行三次Hermite 计算，然后作基线漂移校正 . 具体步骤如下.图2.3 确定插值结点Step1 确定信号长度m.如下确定了几个典型地m 值，如表2-2所示.Step3 计算插值地边界条件.实验是用8组不同心率地心电信号，迭加上不同频率地基线漂移（0.1Hz,0.2Hz10.3Hz)来进行基漂校正.图2.4所示为其中一例，心率为105次/分.迭加三种不同频率地正弦波作为基线漂移.图中1C 为原始信号,2C 为由插值函数计算所得地基漂，3C 为经校正后地心电信号.图2.4 实时基漂校正结果.HR=105次/分.1C 原始ECG ，2C 由插值拟合地基漂，3C 校正后地ECG.而实际临床情况，心率一般均在60次/分以上，基漂频率为0.17Hz 至0.33Hz 之间，所以基漂校正地仿真结果误差都在1.0%以下，可以满足临床要求.ST 段地计算也是令人满意地.硬件电路实现虽然上述插值方法经过改进与简化，计算量已有很大减少，但对小型实时心电监护仪来说，CPU 还可能不能承担.因此又用专用硬件电路实现了上述地插值计算，并且还构成了一个ST 段检测仪.其框图如图2.5所示.其中插值基函数电路是将Hermite 插值基函数)(),(),(),(4321k k k k ϕϕϕϕ，其中（k=12m,m=375).计算并量化后写入EPROM 片.再在乘法控制线地控制下可依次读出.插值条件寄存电路则由由CPU 送入地每段插值结点地边界条件'1'0,10,,y y y y ,它们也可在乘法控制线地控制下依次读出.这样每当由插值基函数电路端口读出函数值时，乘法控制线变回产生含有4个负脉冲地脉冲序列，乘法电路就依次产生4个对应地乘积)()(),(),(4'1312'010k y k y k y k y ϕϕϕϕ和,这四个乘积经累加电路累加后送至输出端口，完成一次基漂拟合值地计算.由此连续运行k=0至m ，即可完成一个周期地基漂校正.这个电路具有高速性能，插值基函数地确定、乘法运算、累加、翻转、技数、清零、重复等操作均是在乘法控制线地控制下同步进行地，有一部分操作室并行进行地，最大限度地减少了运算时间，提高了运算速度，可在12微秒内完成一个点地插值计算，时钟脉冲频率为100MHz.且整个电路地成本也很低.此监测仪对于心率在40次/分以上，基线漂移频率在0.4Hz 以下地ECG 基线漂移能相当好地进行校正.对ST段地检测，在一般情况下也能满足临床要求.图2.5 插值计算硬件电路框图参考文献[1] 文畅平.人民黄河.湖南：邵阳学院，2006.[2] 李庆阳，王能超，易大义.数值分析[M].北京：清华大学出版社，2008.[3] 白峰杉.数值计算引论.北京：高等教育出版社，2004.[4] 李庆阳.计算科学方法基础.北京：清华大学出版社，2006.[5] 冯康等.数值计算方法.北京：国防工业大学，1978.[6] 张雪敏.MA TLAB基础及应用.北京：中国电力出版社，2009.附录A 三次Hermite插值地C程序 1. 流程图2.C程序代码#include<stdio.h>#include<math.h>float f0(float x){return((x-1)*(x-1)*(2*x+1))。

hermite插值计算公式matlab例题Hermite 插值是一种插值方法，它不仅通过已知的函数值来拟合一个多项式，还利用函数在给定点的导数信息。

这使得 Hermite 插值对于需要在给定点处匹配函数值和导数值的问题非常有用。

下面是一个在 MATLAB 中进行 Hermite 插值计算的简单例子。

假设我们有以下数据点：(1, f(1)), (2, f(2)), (1, f'(1)), (2, f'(2))。

我们要使用 Hermite 插值来计算 f(1.5) 的近似值。

首先，我们定义一个 MATLAB 函数来进行 Hermite 插值：function result = hermite_interpolation(x, y, dy, xi)% x: 数据点的 x 坐标% y: 数据点的函数值% dy: 数据点的导数值% xi: 要进行插值的点的 x 坐标n = length(x);% 初始化插值结果result = 0;for i = 1:n% 计算 Lagrange 基函数的权重Li = 1;for j = 1:nif j ~= iLi = Li * (xi - x(j)) / (x(i) - x(j));endend% Hermite 插值基函数的权重Hi = (1 - 2 * (xi - x(i)) * polyval(polyder(poly(x)), x(i))) * Li^2;% 插值结果的累加result = result + y(i) * Hi + dy(i) * Li;endend接下来，我们使用这个函数进行 Hermite 插值：% 给定数据点x = [1, 2];y = [f(1), f(2)];dy = [f_prime(1), f_prime(2)];% 要进行插值的点xi = 1.5;% 进行 Hermite 插值result = hermite_interpolation(x, y, dy, xi);% 显示插值结果disp(['Hermite interpolation at x = ', num2str(xi), ' is ', num2str(result)]);在这里，f 和 f_prime 是你实际的函数和导数函数。