数值计算数值分析实验分段三次埃尔米特hermite插值c程序实现流程图
实验四分段三次埃尔米特插值
(一)实验目的
掌握分段三次埃尔米特插值算法。
(二)实验项目内容
1.写出计算步骤和流程图。
2.对每种算法分别用C或c#程序实现。
3.调试程序。
可用以下数据进行调试。
已知函数y=1/(1+x2)在区间[0,3]上取等距插值节点,求区间[0,3]上的分段三次埃尔米特插值函数,并利用它求出f(1.5)的近似值(0.3075)。
(三)主要仪器设备
微机
(四)实验室名称
公共计算机实验室
(五)实验报告撰写
实验四分段三次埃尔米特插值
实验报告
一、流程图
二、程序代码
#include<stdio.h> #include<math.h>
float f0(float x)
{
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);
}
三、运行结果【截图】。
三次埃尔米特插值
《计算方法》课程设计报告学生姓名:张学阳学号:1009300132陈洋1009300109刘睿1009300122 学院:理学院班级: 数学101题目: 分段线性及三次埃尔米特插值通用程序指导教师:宋云飞职称:讲师朱秀丽讲师尚宝欣讲师2012年12月30日目录目录 (I)一、摘要 (1)二、算法设计 (1)2.1分段线性插值 (1)2.2分段三次埃尔米特插值 (1)2.3功能框图 (1)三、例题计算 (1)四、误差及结果分析 (9)4.1例题误差分析 (1)4.2结点个数对插值结果的影响 (1)五、总结及心得体会 (12)参考文献 (13)源程序 (14)一、摘要分段线性插值与分段定义的线性插值,在相邻插值节点的区间上对应的是同一个线性函数。
由于它们的表现形式不一样从而产生为两种不同的计算方法,相应的误差表现形式也不一样.拉格朗日插值余项利用f(x)的二阶导数,要f(x)的二阶导数存在,对于二阶导数不存在的情况不能估算出它的误差,所以适用范围比较小.现在我们可以利用一阶导数就估算出误差,给计算带来许多的方便。
为了避免高次插值可能出现的大幅度波动现象,在实际应用中通常采用分段低次插值来提高近似程度,比如可用分段线性插值或分段三次埃尔米特插值来逼近已知函数,但它们的总体光滑性较差。
为了克服这一缺点,一种全局化的分段插值方法——三次样条插值成为比较理想的工具。
在代数插值过程中,人们为了获得较好的近似效果,通常情况下是增加插值节点数.由于二次插值比线性插值近似效果好,因此容易错误地认为插值多项式次数越高越好.事实上,随着插值节点的增多,插值多项式不一定收敛到被插值函数.。
通过分段低次插值或样条插值可以得到较好的近似逼近函数,分段低次插值具有公式简单、运算量小、稳定性好、收敛性有保证等优点.随着子区间长度h取得足够小,分段低次插值总能满足所要求的精度.因此分段低次插值应用十分广泛.。
分段线性插值是分段低次插值中常见的方法之一,在本文中对函数在(-5,5)上进行分段线性插值,取不同节点个数n,得到不同分段线性插值函数.并用MATLAB编写分段线性插值函数,最后比较用不同节点数所得插值函数与真实函数的误差,从而得出节点数与插值效果的关系。
分段三次埃尔米特插值
分段三次埃尔米特插值分段线性插值函数的导)(I x h 数是间断的,若在节点k x (k =0,1,…,n )上除已知函数值外还给k f 出导数值k k m f ='(k =0,1,…,n ),这样就可构造一个导数连续的分段插值函数)(I x h ,它满足条件:(1).[][]),(,)(I 11b a C b a C x h ∈代表区间一[]b a ,阶导数连续的函数集合. (2)k k h f x =)(I ,'')(I k k h f x -(k =0,1,…n ). (3))(I x h 在每个小区间上是三次[]1,+k k x x 多项式.由两点三次插值多项式可以知道在)(I x h 区间上表达[]1,+k k x x 式为'21111211211)()(21*)()21()()(I k k k k k k k k k k k k k k k k k k h f x x x x x x x x x x x x x x f x x x x x x x x x ---+--+--+--+--=++--++++)(+'1121)(+++---k k kk k f x x x x x x )(.分段三次埃尔米特插值比分段线性插值效果明显改善,但这种插值要求给出节点上的导数值,所要提供的信息太多,其光滑度也不高(只有一阶导数连续),改进这种以克服其缺点就导致三次样条插值的提出.三次样条插值上面讨论的分段插值函数都有一致收敛性,但光滑性比较差,对于像告诉飞机的机翼形线,船体放样等型值线往往要求有二阶光滑度,即有二阶连续导数.早起工程师制图是,把富有弹性的细长木条(所谓样条)用压铁固定在样点上,在其他地方让它自由的弯曲,然后画下长条的曲线,称为样条曲线.样条曲线实际上有分段三次曲线并接而成,在连接点即样点上要求二阶导数连续,从而数学上加以概括就得到数学样条这一概念.三次样条函数定义 若函数[]b a C x S ,)(2∈,并且在每个小区间上是[]1,+j j x x 三次多项式,其中是给定b x x x a n =<<<= 10节点,则称是节点)(x S 0x ,1x ,…,n x 上的三次样条函数.若在节点上j x 给定函数值)(j j x f y =(j =0,1,…,n ),并且成立j j y x S =)( (j =0,1,…,n ),(1.1) 则称为三次)(x S 样条插值函数.由定义知道要求出)(x S ,在每个小区间上要确定[]1,+j j x x 4个待定系数,一共有个小n 区间,所以应该确定4个参数n .根据在上二)(x S []b a ,阶导数的连续性,在节点j x (j =1,2,…,n -1)处应该满足连续性的条件)0()0(+=-j j x S x S ,)0()0(''+=-j j x S x S (1.2))0()0(''''+=-j j x S x S .一共有3n -3个条件,再加上要满)(x S 足插值条件(1.1),共有4n -2个条件,因此还需要2个条件才能确定)(x S .通常可以在区间[]b a ,端点0x a =,n x b =上各加上一个条件(称为边界条件),可根据实际问题要求给定.常见的有以下3种;(1)已知两端的一阶导数值,即'00')(f x S =,;')(n n f x S =. (1.3)(2)两端的二阶导数已知,即''00'')(f x S =,'''')(n n f x S =, (1.4)其特殊情况为0)()(''0''==n x S x S . (1.5)(3)当)(x f 是以n x -0x 为周期的周期函数时,则要求也是)(x S 周期函数.这时边界条件应满足)0()0(0-=+n x S x S ,)0()0('0'-=+n x S x S , )0()0(''0''-=+n x S x S . (1.6)而此时(1.1)中n y y =0.这样确定的样条函数称)(x S 为周期样条函数. 埃尔米特插值不少实际问题的插值问题不但要求在节点上函数值相等,而且还要求对应的导数值也相等,甚至要求高阶导数也相等,满足这种要求插值的多项式就是埃尔米特(Hermi t e )插值多项式.下面只讨论函数值与导数值个数一样的情况.设在节点上b x x x a n ≤<<<≤ 10,)(i i x f y =,)('j j x f m =(j =0,1,…,n ),要求插值多项式)(x H ,满足条件j j y x H =)(,j j m x H =)('(j =0,1,…,n ). (1.1)这里给出了2n +2个条件,可唯一确定一个次数不超过2n +1的多项式)()(12x H x H n =+,其形式为12121012)(++++++=n n n x a x a a x H .如果根据条件(1.1)来确定2n +2个系数0a ,1a ,…,12+n a ,显然非常复杂,因此,我们依旧采用拉格朗日插值多项式的基函数的方法.先求插值基函数)(x j α及)(x j β(j =0,1,…,n ),一共有2n +2个,每一个基函数都是2n +1次多项式,且满足条件⎪⎪⎭⎪⎪⎬⎫====⎩⎨⎧=≠==).,,1,0,()(,0)(;0)(,,1,,0)(''n k j x x x k j k j x jk k j k j k j jk k j δββαδα (1.2)于是满足条件(1.1)的插值多项式可以写成)()(12x H x H n +=用插值基函数表示的形式[]∑=-+=nj j j j j n x m x y x H 012)()()(βα. (1.3)由条件(1.2)可以知道,有k k n y x H =-)(12,kn m x H =+)('12,(k =0,1,…,n ).下面的问题就是求满足条件(1.2)的基函数以)(x j α及)(x j β.所以,我们可以利用拉格朗日插值基函数)(x l j .令)()()(2x l b ax x a j j +=,由条件(1.2)有1)()()(2=+=j j j j j x l b ax x α,[]0)()(2)()()(''=++=j j j j j j j j j x l b ax x al x l x α,整理得⎩⎨⎧=+=+0)(21'j j j x l a b ax . 解出)(2'j j x l a -=,)(21'j j j x l x b +=.由于)())(()()())(()()(110110n j j j j j j n j j j x x x x x x x x x x x x x x x x x l --------=+-+- ,利用两边取对数再求导数,有∑≠=-=njk k kj j jx x x l 0'1)(,所以有)()1)(21()(20x l x x x x x a j njk k kj j j ∑≠=---=. (1.4)同理,可以得到)()()(2x l x x x j j j -=β. (1.5)同时还证明满足条件(1.1)的插值多项式是唯一的.用反证法,假设及都满)(12x H n +)(12x H n +足条件(1.1),所以有)()()(1212x H x H x n n ++-=ϕ在每个节点上均有二重根,即)(x ϕ有2n +2重根.但是是不高)(x ϕ于2n +1次的多项式,所以0)(≡x ϕ.唯一性得到证明.。
Hermite插值方法
数值分析实验报告五一、实验目的理解Hermite插值方法,掌握Hermite插值算法设计二、实验内容使用vc++编程,实现该方法,即Hermite插值法三、实验步骤#include <iostream.h>double herm(double x0,double x1,double y0,double y1,double h0,double g0,double g1,double x) {double alp0,alp1,bta0,bta1,t;double s;t=h0*h0;alp0=(x-x1)*(x-x1)*(h0+2*(x-x0))/t/h0;alp1=(x-x0)*(x-x0)*(h0-2*(x-x1))/t/h0;bta0=(x-x0)*(x-x1)*(x-x1)/t;bta1=(x-x1)*(x-x0)*(x-x0)/t;s=y0*alp0+y1*alp1+g0*bta0+g1*bta1;return(s);}void main(){int n=7;double p0;double pn; double aa[8],bb[8],s=0;double xx[8]={0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9};double yy[8]={0.4794,0.6442,0.7833,0.8912,0.9636,0.9975,0.9917,0.9463};double g[8];int i;double a[8],c[8],h[8];cout<<"Please input p0 and pn"<<endl;cin>>p0;cin>>pn;for(i=0;i<=n-1;i++){h[i]=xx[i+1]-xx[i];cout<<"h["<<i<<"]="<<h[i]<<endl;}c[0]=1;g[0]=3*(yy[1]-yy[0])/h[0]-p0*h[0]/2;for( i=1;i<=n-1;i++){a[i]=h[i]/(h[i]+h[i-1]);c[i]=1-a[i];}for(i=1;i<n;i++){cout<<"a["<<i<<"]="<<a[i]<<endl;cout<<"c["<<i<<"]="<<c[i]<<endl;}for( i=1;i<=n-1;i++){g[i]=3*(c[i]*(yy[i+1]-yy[i])/h[i]+a[i]*(yy[i]-yy[i-1])/h[i-1]);}a[n]=1;g[n]=3*(yy[n]-yy[n-1])/h[n-1]+pn*h[n-1]/2;for(i=0;i<=n;i++)cout<<"g["<<i<<"]="<<g[i]<<endl;aa[0]=2;bb[0]=c[0]/aa[0];g[0]=g[0]/aa[0];for(i=1;i<=n-1;i++){aa[i]=2-a[i]*bb[i-1];bb[i]=c[i]/aa[i];g[i]=(g[i]-a[i]*g[i-1])/aa[i];}aa[n]=2-a[n]*bb[n-1];g[n]=(g[n]-a[n]*g[n-1])/aa[n];for(i=n-1;i>=0;i--){g[i]=g[i]-bb[i]*g[i+1];}cout<<endl;for(i=0;i<=n;i++)cout<<"g["<<i<<"]="<<g[i]<<endl;double ss;double c0,c1,d0,d1,g0,g1,h1;double x0;cout<<"Please input interpolation point x0:"<<endl;cin>>x0;if(x0>=0.5 && x0<0.7){c0=xx[0];c1=xx[1];d0=yy[0];d1=yy[1];h1=h[0];g0=g[0];g1=g[1];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=0.7 && x0<0.9){c0=xx[1];c1=xx[2];d0=yy[1];d1=yy[2];h1=h[1];g0=g[1];g1=g[2];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=0.9 && x0<=1.1){c0=xx[2];c1=xx[3];d0=yy[2];d1=yy[3];h1=h[2];g0=g[2];g1=g[3];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=1.1 && x0<=1.3){c0=xx[3];c1=xx[4];d0=yy[3];d1=yy[4];h1=h[3];g0=g[3];g1=g[4];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=1.3 && x0<=1.5){c0=xx[4];c1=xx[5];d0=yy[4];d1=yy[5];h1=h[4];g0=g[4];g1=g[5];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=1.5 && x0<=1.7){c0=xx[5];c1=xx[6];d0=yy[5];d1=yy[6];h1=h[5];g0=g[5];g1=g[6];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}else if(x0>=1.7 && x0<=1.9){c0=xx[6];c1=xx[7];d0=yy[6];d1=yy[7];h1=h[6];g0=g[6];g1=g[7];ss=herm(c0,c1,d0,d1,h1,g0,g1,x0);cout<<ss<<endl;}elsecout<<"The data error,please input again!"<<endl;}四、运行结果。
埃尔米特(Hermite)插值
实验二埃尔米特(Hermite)插值一、实验目的:1.掌握埃尔米特插值算法原理;2.使用C语言编程实现埃尔米特插值算法。
二、实验准备:阅读《数值分析》2.4节二、实验要求:某人从甲地开车去乙地,每隔一段时间对行车距离和速率进行一次采样,得到在n+1 个采样时刻点t i 的里程s i和速率v i(i=0, 1, ..., n)。
要求编程构造埃尔米特插值多项式H2n+1(t),满足H2n+1(t i)=s i,H'2n+1(t i)=v i,对所有i=0, 1, ..., n成立,并据此计算m个给定时刻的里程和速率。
函数接口定义:void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] );其中N为采样点个数(注意这个N不是公式中的最大下标n,而是等于n+1),采样时刻点t i、里程s i、速率v i分别通过t、s、v传入;m是需要估算的给定时刻的个数,ht传入给定的时刻点,相应计算出的里程和速率应分别存储在hs和hv中。
裁判程序如下:裁判输入数据:20.0 1.00.0 1.00.0 0.050.0 0.2 0.5 0.8 1.030.0 0.5 1.0100.0 170.0 200.030.0 150.0 0.050.0 0.25 0.5 0.75 1.050.0 1.0 2.0 3.0 4.00.0 60.0 160.0 260.0 300.05.0 70.0 100.0 120.0 20.0100.5 1.0 1.5 2.0 2.5 3.0 3.5 3.8 3.95 4.0标准输出数据:0.0000 0.1040 0.5000 0.8960 1.00000.0000 0.9600 1.5000 0.9600 0.0000100.0000 127.9297 170.0000 195.9766 200.000030.0000 165.4688 150.0000 52.9688 0.000030.2222 60.0000 105.9303 160.0000 206.3438 260.0000 307.9764 305.7687 299.9796 300.000062.6024 70.0000 109.0488 100.0000 92.9745 120.0000 41.2374 -44.8421 -16.2783 20.0000#include<stdio.h>#define MAXN 5 /* 最大采样点个数 */#define MAXM 10 /* 最大估算点个数 */void Hermite_Interpolation( int N, double t[], double s[], double v[], int m, double ht[], double hs[], double hv[] ){double l[10],p[10],h1[10],h2[10],x,ll[10],pp[10];int kk;for(kk=0;kk<m;kk++){x=ht[kk];hs[kk]=0;hv[kk]=0;int i;for(i=0;i<N;i++){l[i]=1;ll[i]=1;int j;for(j=0;j<N;j++){if(i!=j){l[i]=l[i]*(x-t[j])/(t[i]-t[j]);}}p[i]=0;pp[i]=0;int k;for(k=0;k<N;k++){if(i!=k){p[i]=p[i]+l[i]/(x-t[k]);pp[i]=pp[i]+ll[i]/(t[i]-t[k]);}}h1[i]=(1-2*pp[i]*(x-t[i]))*l[i]*l[i];h2[i]=(x-t[i])*l[i]*l[i];hs[kk]=hs[kk]+s[i]*h1[i]+v[i]*h2[i];int kkk;for(kkk=0;kkk<N;kkk++){if(x==t[kkk])break;}if(x==t[kkk])hv[kk]=v[kkk];elsehv[kk]=hv[kk]+s[i]*(2*p[i]*l[i]-4*l[i]*p[i]*(x-t[i])*pp[i]-2*pp[i]*l[ i]*l[i])+v[i]*(l[i]*l[i]+2*l[i]*p[i]*(x-t[i]));}}}int main(){int N, m;double t[MAXN], s[MAXN], v[MAXN]; /* 用于构造的数据 */double ht[MAXM], hs[MAXM], hv[MAXM]; /* 用估算的数据 */int i;while ( scanf("%d", &N) != EOF ) {for ( i=0; i<N; i++ )scanf("%lf", &t[i]);for ( i=0; i<N; i++ )scanf("%lf", &s[i]);for ( i=0; i<N; i++ )scanf("%lf", &v[i]);scanf("%d", &m);for ( i=0; i<m; i++ )scanf("%lf", &ht[i]);Hermite_Interpolation( N, t, s, v, m, ht, hs, hv );for ( i=0; i<m; i++ )printf("%.4lf ", hs[i]);printf("\n");for ( i=0; i<m; i++ )printf("%.4lf ", hv[i]);printf("\n\n");}return 0; }。
埃尔米特(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; }。
分段线性及三次埃尔米特插值通用程序
结 束
计算方法课程设计
分段线性及三次埃尔米特插值通用程序
-----------数学101 张学阳 陈 洋 刘 睿
分段线性插值在每个小区间上可表示为
I h ( x) x xk 1 x xk f k f k 1 , xk x xk 1 , k 0,1,, n 1. xk xk 1 xk 1 xk
M = -5:0.01:5; output = zeros(1,length(M)); n = 1; for i=2:N+1 for x=-5:0.01:5 if x<xx(i) && x>=xx(i-1) lx(1)=ff(i-1)*(x-xx(i))/(xx(i-1)-xx(i)); lx(2)=ff(i)*(x-xx(i-1))/(xx(i)-xx(i-1)); output(n) = lx(1)+lx(2); n = n+1; end end end
end
end end ezplot(f,[-5,5]) hold on A =-5:0.01:5; plot(A,output,'r'); end if B~=1&&B~=2 disp('输入有误,请重新输入') end c=input('是否继续(1为继续):'); close end
M = -5:0.01:5; output = zeros(1,length(M)); n = 1; for i=2:N+1 for x=-5:0.01:5 if x<xx(i) && x>=xx(i-1) lx(1)=ff(i-1)*(x-xx(i))^2/(xx(i-1)-xx(i))^2*(1+2*(x-xx(i1))/(xx(i)-xx(i-1))); lx(2)=ff(i)*(x-xx(i-1))^2/(xx(i)-xx(i-1))^2*(1+2*(xxx(i))/(xx(i-1)-xx(i))); lx(3)=ff1(i-1)*(x-xx(i))^2/(xx(i-1)-xx(i))^2*(x-xx(i-1)); lx(4)=ff1(i)*(x-xx(i-1))^2/(xx(i)-xx(i-1))^2*(x-xx(i)); output(n) = lx(1)+lx(2)+lx(3)+lx(4); n = n+1;
数值分析实验六(分段三次Hermite插值)
数值分析实验六(分段三次Hermite插值)《数值分析》实验报告实验编号:实验六课题名称:分段三次Hermite插值一、算法介绍给定的函数为f(x)=1/(25*x*x+1),将给定区间分成10分,得到11个节点:x[0],x[1],...,x[10],构造插值函数的基函数。
当x在(x[0],x[1])区间上时,H[0] = (x-x[0])*[((x-x[1])/(x[0]-x[1]))^2]。
其余的区间为H[0]=0。
h[0]= [1+2*(x-x[0])/(x[1]-x[0])]*[((x-x[1])/(x[0]-x[1]))^2]。
当x在[x[i-1],x[i]] (i=1,2,3, (9)区间上时,H[i]=(x-x[i])*[((x-x[i-1])/(x[i]-x[i-1]))^2],h[i]=[1+2*(x-x[i])/(x[i-1]-x[i])]*[((x-x[i-1])/(x[i]-x[i-1]))^2)。
当x在(x[i],x[i+1]](i=1,2,3,…,10)区间上其余的区间为H[i]=(x-x[i])[((x-x[i+1])/(x[i]-x[i+1]))^2],h[i]=[1+2*(x-x[i])/(x[i+1]-x[i])]*[((x-x[i+1 ])/(x[i]-x[i+1]))^2]。
其余区间上均为H[i]=0,h[i]=0(i=1,2,…,10)。
当x在(x[9],x[10])区间上时,H[10] = (x-x[9])(((x-x[10])/(x[9]-x[10]))^2).其余的区间为H[10]=0.h[10]= (1+2*((x-x[9])/(x[10]-x[9])))(((x-x[10])/(x[9]-x[10]))^2).其余区间h[10]=0。
构造函数H(x) =∑(y[i]*h[i]+y'[i]*H[i],(i=0,1,…,10)。
二、程序代码// testV iew.cpp : implementation of the CT estV iew class//#include "stdafx.h"#include "test.h"#include "testDoc.h"#include "testView.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////////////////////////////// //////////// CTestV iewIMPLEMENT_DYNCREA TE(CTestView, CView)BEGIN_MESSAGE_MAP(CTestView, CView)//{{AFX_MSG_MAP(CTestView)// NOTE - the ClassWizard will add and remove mapping macros here.// DO NOT EDIT what you see in these blocks of generated code!//}}AFX_MSG_MAP// Standard printing commandsON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_DIRECT, CV iew::OnFilePrint) ON_COMMAND(ID_FILE_PRINT_PREVIEW,CView::OnFilePrintPreview)END_MESSAGE_MAP()/////////////////////////////////////////////////////////////////// //////////// CTestV iew construction/destructionCTestView::CTestV iew(){// TODO: add construction code here}CTestView::~CT estView(){}BOOL CTestView::PreCreateWindow(CREA TESTRUCT& cs){// TODO: Modify the Window class or styles here by modifying // the CREA TESTRUCT csreturn CV iew::PreCreateWindow(cs);}/////////////////////////////////////////////////////////////////// //////////// CTestV iew drawingvoid CTestView::OnDraw(CDC* pDC){CTestDoc* pDoc = GetDocument();ASSERT_V ALID(pDoc);// TODO: add draw code for native data hereint i,j,k;double x,y,p_x,p_y,l,xx[100],f[100],F[100],sum,p_sum;CPen MyPen,*OldPen;pDC->SetViewportOrg(400,400); //定义坐标原点for(i=-500;i<500;i++){pDC->SetPixel(0,i,RGB(0,0,0));pDC->SetPixel(i,0,RGB(0,0,0)); //画出坐标}pDC->TextOut(-210,5,"-1");pDC->TextOut(196,5,"1");//原函数MyPen.CreatePen(PS_SOLID,1,RGB(255,0,0));//定义画笔颜色OldPen=pDC->SelectObject(&MyPen);x=-1.0,y=1/(1+25*x*x);p_x=x*200;p_y=-y*200;pDC->MoveTo(p_x,p_y);for (x=-1.0;x<=1.0;x+=0.0001){y=1/(1+25*x*x);p_x=x*200;p_y=-y*200;pDC->LineT o(p_x,p_y);}pDC->SelectObject(OldPen);MyPen.DeleteObject();//分段三次Hermite插值MyPen.CreatePen(PS_SOLID,1,RGB(0,0,0)); OldPen=pDC->SelectObject(&MyPen); x=-1.0,y=1.0/(1+25*x*x);p_x=x*200;p_y=-y*200;pDC->MoveTo(p_x,p_y);x=-1.0;for(i=0;i<=10;i++){f[i]=1/(1+25*x*x);xx[i]=x;F[i]=-50*x/(1+25*x*x)/(1+25*x*x); //导数x+=0.2;}x=-1.0;for(j=0;j<=1000;j++){sum=0;for(i=0;i<=10;i++){if(x==xx[i]){sum=f[i];p_x=x*200,p_y=-sum*200;pDC->LineT o(p_x,p_y);break;}if(xxx[i]){y=(1+2*(x-xx[i])/(xx[i+1]-xx[i]))*(x-xx[i+1])*(x-xx[i+1])/(xx[i]-xx[i+1])/(xx[i]-xx[i+1]);sum+=f[i]*y;y=(1+2*(x-xx[i+1])/(xx[i]-xx[i+1]))*(x-xx[i])*(x-xx[i])/(xx[i+1]-xx[i])/(xx[i+1]-xx[i]);sum+=f[i+1]*y;y=(x-xx[i])*(x-xx[i+1])*(x-xx[i+1])/(xx[i]-xx[i+1])/(xx[i]-xx[i+1]);sum+=F[i]*y;y=(x-xx[i+1])*(x-xx[i])*(x-xx[i])/(xx[i+1]-xx[i])/(xx[i+1]-xx[i]);sum+=F[i+1]*y;p_x=x*200;p_y=-sum*200;pDC->LineT o(p_x,p_y);break;}}x+=0.002;}pDC->SelectObject(OldPen);MyPen.DeleteObject();/////////////////////////////////////////////////////////////////// //////////// CTestV iew printingBOOL CTestView::OnPreparePrinting(CPrintInfo* pInfo){// default preparationreturn DoPreparePrinting(pInfo);}void CTestView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add extra initialization before printing}void CTestView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add cleanup after printing}/////////////////////////////////////////////////////////////////// //////////// CTestV iew diagnostics#ifdef _DEBUGvoid CTestView::AssertV alid() const{CView::AssertV alid();}void CTestView::Dump(CDumpContext& dc) const{CView::Dump(dc);CTestDoc* CT estV iew::GetDocument() // non-debug version is inline{ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CT estD oc)));return (CT estDoc*)m_pDocument;}#endif //_DEBUG/////////////////////////////////////////////////////////////////// //////////// CTestV iew message handlers三、运算结果截屏红色的曲线为原函数图像,黑色曲线为分段三次Hermite插值曲线四、算法分析上述图像中黑色的曲线为分段分段三次Hermite插值多项式所对应的图像,由图像可看出黑色的分段三次Hermited插值函数图像和拉格朗日、分段线性插值相比与红色被逼近函数的重合度最好,说明分段三次Hermite插值在函数的各节点两边插值函数的导数是相等的,保证了在各节点的平滑性,且不会出现Runge现象。
埃尔米特插值
2013-2014(1)专业课程实践论文题目:埃尔米特插值一、算法理论1、埃尔米特插值多项式设已知函数()y f x =在节点n x x x <<< 10上的函数值),,1,0)((n i x f y i i ==以及一切导数值()(0,1,)i i y f x i n ''==,要求一个插值多项式()H x ,使其满足()i i H x y =,()i i H x y ''= ,n i ,,1,0 = (1)显然,由条件(1)可以确定一个次数不高于21n +的代数多项式)(12x H n +,曲线)(12x H y n +=与()y f x =在节点处不仅重合而且有公共切线。
我们采用拉格朗日插值基函数的方法。
先求插值基函数),,1,0)((),(n j x x j j =βα,共22n +个基函数,每一个基函数都是一个21n +次多项式,且满足条件(),()0;()0,(),0,1,,.j jk j k j k j k jk x x x x j n αδαββδ'==='== (2)这里⎩⎨⎧=≠=.,1,,0k j k j δ(3)于是满足条件(1)的插值多项式 可写成用插值函数表示的形式210()[()()]nn j j j j H x y x y x αβ+='=+∑(4)由条件(2),显然有2121(),()(0,1,,).n i i n i j H x y H x y i n ++''===下面的问题就是要求满足条件(2)的)(x j α与).(x j β为此,可利用拉格朗日插值基函数)(x l j ,由条件(2)有n 个二重零点),,,1,0(j k n k x k ≠= ,于是可令).()()(2x l b ax x j j +=α 由条件(2)有2()()()1,()()[()2()()]0.j j j j j j j j j j j j j x ax b l x x l x al x ax b l x αα=+=''=++=解出2(),12().j j j j j a l x b x l x ''=-=+由于,)(11111100nj nj j j j j j j j j x x x x x x x x x x x x x x x x x x x x x l ----------=++--故01(),nj j k j kk jl x x x =≠'=-∑于是).(]1)(21[)(20x l x x x x x j njk k kj j j ∑≠=---=α (5)同理,可得).()()(2x l x x x j j j -=β (6) 将(5)式、(6)式 代入式(4)便得到埃尔米特插值多项式22210001()[12()]()()()nnnn i i j j j j j k j j kk jH x y x x l x y x x l x x x +===≠'=--+--∑∑∑(7)满足条件(1)的埃尔米特插值多项式是唯一的。
5.4 三次Hermite插值
由于0 ( x0 ) 0 ( x1 ) 0( x1 ) 0, 故 0 ( x) 含有 ( x x0 )( x x1 ) 因子。可设
0 ( x) c( x x0 )( x x1 )2
2
其中c为待定系数。
1 . 由 0( x0 ) 1, 可得 c 2 ( x0 x1 )
称 H3 ( x) 为三次Hermite插值多项式,上述条件称 为此问题的插值条件。
采用基函数的方法来构造 H3 ( x) 。
将 H3 ( x) 表示为:
H3 ( x) y00 ( x) y11( x) m00 ( x) m11( x)
其中 0 ( x),1 ( x), 0 ( x), 1 ( x)为插值基函数,且均为次数 不超过3的多项式。为满足插值条件,它们应满足
并给出余项公式。
( x1 ) m1, (i 0,1,2), H3 ( xi ) yi , H3 H3 ( x) y00 ( x) y11( x) y22 ( x) m11( x)
函数值
导数值
x0
x1
0
1
x2
0 0
1
x1
0 0 0
1
0 ( x)
1( x)
2 ( x )
作为多项式插值,三次已是较高的次数,次数再高就有 可能发生Runge现象,因此,对有n+1节点的插值问题, 我们可以使用分段两点三次Hermite插值
例1 对给定数表
x
y
x0 y0
x1 y1 m1
x2 y2
y
求一个三次Hermite插值多项式 H3 ( x) 满足条件
( x1 ) m1, (i 0,1,2), H3 ( xi ) yi , H3
三次Hermite插值
04 实例分析
CHAPTER
实例一:已知数据点的插值
总结词
利用已知数据点进行插值,可三次Hermite插值方法,利用已知的数据点来估计未知点的值。这 种方法能够更好地处理数据点的变化,并提高插值的精度。
CHAPTER
插值多项式的构造
定义
Hermite插值法是一种通过已知的离散数据点来构造一个多 项式,使其能够准确地经过这些数据点,并尽可能地平滑地 连接这些点的方法。
构造方法
Hermite插值多项式由两个部分组成,一个是线性函数,另 一个是二次函数。线性函数部分用于确保插值多项式能够准 确地经过数据点,而二次函数部分则用于保证插值多项式的 平滑性。
实例二:未知数据点的插值
总结词
在未知数据点的情况下,可以通过三次 Hermite插值方法,预测并估计未知点的值。
详细描述
在数据点未知的情况下,可以利用三次 Hermite插值方法,根据已知的数据点来预 测和估计未知点的值。这种方法能够为后续 的数据分析和处理提供重要的参考依据。
实例三:复杂函数的插值
三次Hermite插值能够提供高精度的插值结果,特别是在处理
复杂函数时。
稳定性好
02
该方法在处理大数据集时表现出良好的稳定性,不易受到噪声
和异常值的影响。
易于实现
03
三次Hermite插值的算法相对简单,易于在计算机上实现和优
化。
三次Hermite插值的局限性
对初始数据敏感
三次Hermite插值的结果对初始数据的选择 较为敏感,不同的初始数据可能导致不同的 插值结果。
