邰艳艳靳由尹(3).doc

数值分析中插值与拟合方法的探讨邰艳艳 靳由尹（2007051322，2007051431）指导教师姓名：郭德龙职 称：讲师单 位：数学系专 业 名 称：信息与计算科学论文提交日期：2011年4月论文答辩日期：2011年4月学位授予单位：黔南民族师范学院答辩委员会主席:论 文 评 阅 人:2011年 月 日33数值分析中插值与拟合方法的探讨邰艳艳，靳由尹（2007051322，2007051331）（黔南民族师范学院数学系信息与计算科学专业 贵州 都匀 558000）摘 要 通过对数值分析中插值与拟合方法的探讨 ,讨论了插值函数的基本概及线性插值和多项式插值存在唯一性 ,主要介绍了基于基函数的拉格朗日插值、基于均差的牛顿插值和基于导数埃尔米特插值.曲线拟合及基于最小二乘拟合的多项式插值.对给出的一组数据(xi,yi)(i=0,1,n),根据不同的原则我们可以利用插值法和曲线拟合法分别来构造近似函数表达式。并通过例子来对以上三种插值方法和基于最小二乘拟合的多项式插值方法进行分析和比较,讨论了各种方法的优点、缺点，从而得出了什么情况下选用何种方法更好。关键词 拉格朗日插值;牛顿插值;埃尔米特插值;曲线拟合;最小二乘法The study of function fitting and interpolation（Mathematics QianNan Normal College for Nationalities, Duyun，Guizhou，558000） Tai yan-yan ，Jin you-yin Abstract Through numerical analysis of the interpolation and fitting methods are discussed,This paper discusses the basic concept of interpolation function and the uniqueness existed in the linear interpolation and polynomials interpolation Mainly introduced based on base function based on the Lagrange interpolation deiation Newton interpolation and based on cubic-hermite interpolating derivative. Curve fitting and based on least squares fitting polynomial interpolation on all three. And the interpolation method and based on the least squares fitting polynomial interpolation method are analyzed and compared.Given a set of data on.According to different principle, we can use the interpolation method and curve-fitting method to construct the approximate function expression respectively, and through examples to to the above three interpolation method and least squares fitting polynomial interpolation method are analyzed and compareddiscusses the advantages and the disadvantages of the different methods to work quantity and its graphics and then work out which method is the best way under the certain circumstance.Key words Lagrange interpolation ; Newton interpolation ; Hermite interpolatiothe curve fitting ; least square fit1 引 言 函数常被用来描述客观事物变化的内在规律（数量关系）.但在生产和科研实践中遇到的大量函数,却是复杂函数.对于实际中的这些复杂函数,我们希望能构造一个能反映函数本身的特性,又便于计算的简单函数,近似代替原来的函数.解决上述问题的方法有两类:一类是对于一组离散,选定一个便于计算的函数形式,如多项式函数、分段性函数、有理函数、三角函数等,要求简单函数满足.由此确定函作为的近似函数,这就是插值方法.另一类方法在选定近似函数的形式时,不要近似函数必须满足,而只要在某种意义下(最小二乘法原理),使近似函数在这些点上的总偏差量最小,这种方法称为曲线拟合.2 插值问题与插值多项式定义1 设为区间上函数，为上互不相同的点，为给定的某一函数类，若上有函数，满足 . (1)则称为关于节点在上的插值函数，称点为插值节点；称为插值型值点，简称型值点或插值点;称为被插函数.定义2 已知函数在区间上的个点的值，即已知，寻求一个解析形式的函数，使之满足 (2)则称为插值结点,为被插值函数，为插值函数，称条件为插值条件，若为次数不超过的多项式，即， . (3)其中为实数，则称为插值多项式.定理1在个相异结点满足插值条件而次数不高于的多项式是唯一的.2.1 拉格朗日插值多项式给定，构造次数不超过的拉格朗日插值多项式 . (4)称为关于的次拉格朗日插值多项式，它满足 . (5)其中称为为结点的次插值函数，它满足 i=0,12,3.n (6). (7)设是上关于的次插值多项式，在上有阶连续导数，在上存在，则其余项为 . (8)例1：根据下表所列的数据点求出其拉格朗日插值多项式，并计算当时的值。表40.84150.9320 0.9738 0.5985 -0.7568解：编写拉格朗日多项式插值方法函数（H:毕业设计程序 Lagrange.m） x=1 1.2 1.8 2.5 4; y=0.8415 0.9320 0.9738 0.5985 -0.7568; f = Lagrange(x,y) f =1.05427*t-.145485e-1*t2-.204917*t3+.328112e-1*t4-.261189e-1 f = Lagrange(x,y,1.6)f =0.9992结果分析：表格中的数据点是按给出的，从插值函数在的值看出，插值函数的精度是比较高的。例2：已知，使用拉格朗日多项式插值法计算，并给出插值多项式。(1) 编写拉格朗日多项式插值方法函数（H:毕业设计程序 lag.m）(2)编写主函数（H:毕业设计程序 lag_main.m）(3)运行程序。 clear lag_main角度du = -40 47 53 79 174插值结果yt = 0.5387 0.6820 0.6017 0.1915 -1.1121cos函数值yreal = 0.7660 0.6820 0.6018 0.1908 -0.9945插值与函数值误差dy = -0.2273 -0.0000 -0.0001 0.0006 -0.1176yt = .136*x3-.66*x2+.9e-1*x+.984 图1结果分析：该程序在输入时需要插值计算的点可以是一序列点，这样程序运行一次就可以同时把这些结果计算出来，避免重复计算，又方便操作实现。程序直接给出插值多项式的表达式：yt = .136*x3-.66*x2+.9e-1*x+.984计算插值结果与函数结果做了比较，对应误差：-0.2273，-0.0000，-0.0001，0.0006，-0.1176。2.2 牛顿插值定义3 零阶均差 一阶均差 . (9)二阶均差. (10)2阶均差是1阶均差的均差，可递推阶均差，得. (11)2.2.1 均差（差商）的性质（）阶均差与函数值的关系为. (12)（）均差关于所含结点是对称的，若为的任意排列，则 (13)即均差值与结点次序无关.2.2.2 牛顿插值多项式给定，次数不超过的牛顿插值多项式为 (14). 表2 (均差表) 2.2.3 插值余项 . (15)由插值多项式的唯一性知，因此，牛顿插值与拉格朗日插值有相同的余项表达式，即 (16) 例3：（1）已知零阶Bessel函数在若干点处得函数值如下： 表31.02.20.76519770.62008600.45540220.28181860.1103623 牛顿插值法计算在1.5处的近似值。（2）5阶勒让德多项式为在-1,1上任意取9个点如=-1，-0.3，0.4，-0.7，0，-0.4，0.8，0.7，1不在乎取点顺序。a） 根据5阶勒让德多项式，算出这9个点的值。b） 根据这9个点作牛顿插值，画出插值多项式图像，并与勒让德多项式图像比较。c） 依据牛顿插值，计算出=0.24，-0.46，0.38时的近似结果，并与真实结果比较。解：编写牛顿插值方法函数（H:毕业设计程序niudun.m）(2)编写主函数（H:毕业设计程序 niudun_main.m）(3)运行程序。 niudun_main运行结果： yt = 0.5118 t = 0.2400 -0.4600 0.8300legen5 = 0.3353 -0.1730 -0.3449chazhi = 0.3353 -0.1730 -0.3449wucha = 1.0e-013 * 0.0061 -0.0028 -0.1260 niudun_main yt = 0.5118t = 0.2400 -0.4600 0.8300legen5 = 0.3353 -0.1730 -0.3449chazhi = 0.3353 -0.1730 -0.3449wucha = 1.0e-013 * 0.0061 -0.0028 -0.1260 图2结果分析：（1）零阶Bessel函数在自变量为x=1.5时，其真值为0.512857，可以看出计算精度还是高。（2）关于勒让德多项式插值问题，对于我们给定的一组要求计算数据，从输出的最后一行，可以看到误差在10的-14到-16次方，这个精度是非常高的。2.2.4 等距结点的牛顿插值 若插值结点为等距结点，即，称为步长，表示在上的值，则有等距结点的牛顿插值公式.定义4 令分别称为在点的一阶向前差分和一阶向后差分。由此可递推阶向前差分和阶向后差分为. (17)并规定零阶差分为 (18)均差与差分有以下关系，即. (19). (20)d） 表4 (差分表)2.2.5 牛顿前插、后插插值公式及其余项 牛顿前插公式为 . (21) 其余项为 . (22) 牛顿后插公式为 .(23) 其余项为 (24) . 例4：根据下表所列的数据点求出其差分形式的牛顿插值多项式，并计算当时的值。 表51.80.84150.9320 0.9854 0.99960.9738 解：编写等距节点插值方法函数（H:毕业设计程序 Newtonforward.m与Newtonback.m） x=1:0.2:1.8; y= 0.8415 0.9320 0.9854 0.9996 0.9738; f = Newtonforward(x,y) f = .841500+.108025*t-.169042e-1*t2-.675000e-3*t3+.541667e-4*t4 f = Newtonback(x,y) f = 973800-.457417e-1*t-.198042e-1*t2+.191667e-3*t3+.541667e-4*t4 f = Newtonforward(x,y,1.55) f = 0.9998 f = Newtonback(x,y,1.55) f = 0.9998 结果分析：表格中的数据点是按给出的，从插值函数在的值看出，牛顿差分插值函数是很准确的。2.3 埃尔米特插值多项式 插值多项式除了满足插值条件外，还要求与被插函数在上面等式共有个条件可唯一确定次数不超过次的多项式，称之为埃米尔特插值多项式，它用插值基函数可表示为 (25)其中，和是插值基函数. 2.3.1 插值基函数 和是满足下列条件 (26)的次多项式.容易求得 (27)其中，是拉格朗日插值基函数.若在插值区间内存在阶导数，则次艾尔米特插值多项式的余项为 (28)其中 (29) 特别地，当时，结点为满足条件 (30)的艾尔米特插值多项式为 (31) (32) . (33) (34)误差估计. (35) 例5：（1）人为的制造一个函数。任意取点列1 、3、 -8、 6 、-4（不需要按照大小顺序排列），计算这些点的函数值，在计算这些点上的导数值。有了这三组值以后，利用Hermite插值方法对以下一些点（或点列）的情况进行插值计算。a） 矩阵中的每个元素。b） 向量X=（-5.3 2.4 -4.2 -1.8 3.4）的每个元素。（2）人为的制造一个函数，计算出这个函数在以下点列中的函数值与导数值。X=（-3，-2，-1，0，1,2,3）利用上面的三组数据插值计算以下点列的函数近似值。P=（-5.1，,3.2，,1.4,1.6,5.3）算出这些点的真函数值，比较不同点的误差情况。解：（1）编写hermite插值的方法函数（H:毕业设计程序herm.m） （2）编写主函数（H:毕业设计程序 herm_main.m）（3）运行程序： format short herm_main回车得到： z1 = x2 z2 = 81 100 625 5625 z3 = 28.0900 5.7600 17.6400 3.2400 11.5600 chazhi = 25.5247 0.5978 -1.9315 2.5589 -26.2206zhen = 24.0804 0.5978 -1.9315 2.5589 -23.3784wucha = 1.4443 0.0000 0.0000 -0.0000 -2.8422图3 图4结果分析：从输出结果来看，结果是让我们非常满意的，这是因为我们对插值的对象是多项式，如果是其他函数，甚至如果不是函数的离散点列可能就没有这么好的效果了。 这也可以看出Hermite插值结果是很不错的。图（3）给出了第一题中插值的函数图象，其实就是二次抛物线。图形中的圆点是插值结果，二次曲线是原来的函数图形。图形中也直接反映了这些点都分布在曲线上。图（4）中，曲线是原函数图象，而*则是离散的插值点，可以看到咋爱区间内符合精度非常高，肉眼很难分辨，而在区间外误差会迅速扩大，甚至很快失效。 3 曲线拟合的最小二乘3.1 最小二乘原理 设已知某物理过程的一组观测数据 ， . （36） 要求在某特定函数类寻求一个函数作为的近似函数，使得二者在 上的误差或称残差 ， （37） 按某种度量标准为最小，这就是拟合问题.要求残差按某种度量标准为最小，即要求由残差构成的残差向量的某种范数为最小，要求，或即 (38)为最小，这本来都是很自然的，可是计算不太方便.以通常要求： (39)或者 （40）为最小.种要求误差平方和最小的拟合称为曲线拟合的最小二乘法.就是说，最小二乘法提供了一种数学方法，利用这种方法可以对实验数据实现在最小平方误差意义下的最好拟合.用最小二乘法求拟合曲线时，必须选择函数类，确定拟合函数的形式，这与所讨论问题的专业知识和经验有关.常，其中，是一组线性无关且已给定的函数，表示组成的函数空间，表示为 （41）此时为线性拟合模型，否则当关于某个或某些参数非线性时，称之为非线性模型.下面给出求解方法.确定出拟合参数，就可得到拟合函数.对于给定的数据，要在给定的函数空间中找一个函数 （42)使满足 (43) 这种求拟合函数的方法称为曲线拟合的最小二乘法.称为最小二乘问题的最小二乘解.令性能指标函数 （44）要使达到极小，由多元函数取极值的必要条件 （45）可得方程组 (46) (47)引入记号 （48）则所得方程组表示成， （49）这个方程组成为正则（或正规）方程组或法方程组，写成矩阵形式为 （50）这是一个系数矩阵为对称性方程组.线性无关时，有唯一解 ， (51)并且相应的拟合函数 （52）就是满足残差平方和为最小的最小二乘解.在实际问题中得到的观测数据并非是等精度、等重要性的.衡量数据的精度和重要性，常常对数据进行“权”处理，对精度好、重要的数据给予较大的权，否则给予小的权，这就是加权最小二乘法.用加权最小二乘法进行拟合是对于观测数据，要求在某函数类中寻求一个函数，使 （53）为最小.中为一组正数，反映数据特性的权，此时正则方程组仍如式(51)，即 只是其中 （54）由以上讨论可知： 对于给定数据，在函数空间中存在唯一函数 (55)使残差平方和为最小. 用最小二乘解得系数可通过解正则方程（51）求得. 用最小二乘解来拟合数据，的平方误差为 (56)例6 ：求数据 表6x-3-2-10123y4230-1-2-5的最小二乘拟合，拟合函数为解：（1）编写主函数（H:毕业设计程序zx_nh）（2）编写拟合的基函数（H:毕业设计程序zx_nh_f）（3）运行程序：clearclczx_nh回车得到： zx_nhc = 0.6667 -1.3929 -0.1310 图5结果分析：上图给出了原离散数据与拟合函数的图象，数据拟合与插值的区别是拟合不需要函数经过原来数据的每一点，而是要求偏差平方和最小。从图形可以看出，在二次多项式拟合的情况下这个表达式在有限的精度内是最优的。 关于基函数的选择，一般是根据所给的数据的曲线图象来判断基函数的特征，有时候则是根据实际经验判断要选择什么样的函数。 与插值方法类似，在不知道函数形态的情况下，有时候如果选择正交多项式作为拟合的基函数效果往往会很不错。 结束语本课题介绍了复杂函数的两种近似表示形式，多项式插值与拟合及其求法.多项式插值是根据插值条件多项式,多项式拟合是根据最小二乘原理,即“偏差平方和最小”求多项式.拉格朗日适合求固定节点的函数值.牛顿插值可以在计算过程中,根据精度要求逐步增加节点,且计算量小,埃尔米特插值适合导数已知的情况,所有插值多项式的次数都不能太高,否则会出现龙格现象.用最小二乘原理进行多项式拟合,适用于根据大量观测数据,寻找变量之间的近似函数关系式.参考文献1 王新和,程世洲,曲线拟合的最小二乘法J，疆职业报,2004,12(2):84-86.2 李桂成，计算方法M,北京,电子工业出版社,2005.10.3 杨泮池,计算方法M,西安,西安交通大学出版社,2005.11.4 马东升,雷勇军,数值计算方法M,2版,北京,机械工业出版社,2009.5 张韵华,奚梅成,陈效群,数值计算方法与算法M,2版,北京,科学社,2006.6 Berden R L, Faires J D, Reynolds A C. Numerical Analysis. Apline Press. 19847 李庆扬 王能超 易大义. 数值分析（第四版）M. 北京清华大学出版社 ,2008. 8 吴才斌. 插值法及其应用J. 湖北大学成人教育学院学报 1999年10月17（5）：9 姜琴，周天宏. 常见的插值法及其应用J. 郧阳师范高等专科学校学报2006年6月第2期10 宋叶志等. MATLAB数值分析与应用M. 北京机械工业出版社,2002.11 维普信息资源系统V6.33, 9/index.asp致谢 经过两个多月的努力，毕业设计终于完成 在整个设计过程中，出现过很多的难题，但都在老师和同学的帮助下顺利解决了，能够顺利完成本次设计，我们首先要各位同学，各位亲朋好友，各位老师，尤其指导老师郭老师，表示忠心的感谢！从确定题目、拟订提纲、完成初步设计到最终的设计完成，我们都得到了老师的精心指导，同学的热心帮助使我们很快掌握了设计的步骤和方法，以及论文的写作方法，并在较短的时间内完成了设计。在不断的学习过程中我体会到：做设计是一个不断学习的过程，从最初刚写论文时对知识问题的模糊认识到最后能够对该问题有深刻的认识，我体会到实践对于学习的重要性，以前只是明白理论，没有经过实践考察，对知识的理解不够明确，通过这次设计，真正做到理论实践相结合。 附录例1：（H:毕业设计程序 Lagrange.m）function f = Lagrange(x,y,x0)syms t;if(length(x) = length(y) n = length(x); else disp(x和y的维数不相等); return;end f = 0.0;for(i = 1:n) l = y(i); for(j = 1:i-1) l = l*(t-x(j)/(x(i)-x(j); end; for(j = i+1:n) l = l*(t-x(j)/(x(i)-x(j); end; f = f + l; simplify(f); if(i=n) if(nargin = 3) f = subs(f,t,x0); else f = collect(f); f = vpa(f,6); end endend例2：（H:毕业设计程序 lag.m）function s=lag(x,y,t)syms p;n=length(x); s=0;for(k=1:n) la=y(k); for(j=1:k-1) la=la*(p-x(j)/(x(k)-x(j); end; for(j=k+1:n) la=la*(p-x(j)/(x(k)-x(j); end; s=s+la; simplify(s); end if(nargin=2) s=subs(s,p,x); s=collect(s); s=vpa(s,4); else m=length(t); for i=1:m temp(i)=subs(s,p,t(i); end s=temp;end(2)（H:毕业设计程序 lag_main.m）x=pi/4,pi/6,pi/3,pi/2;y=cos(pi/4),cos(pi/6),cos(pi/3),cos(pi/2);t=-40*pi/180,47*pi/180,53*pi/180,79*pi/180,174*pi/180;disp(角度)du=-40 47 53 79 174disp(插值结果)yt=lag(x,y,t)disp(cos函数值)yreal=cos(-40*pi/180)cos(47*pi/180)cos(53*pi/180)cos(79*pi/180)cos(174*pi/180)disp(插值与函数值误差)dy=yt-yrealyt=lag(x,y)ezplot(yt,-pi/4,pi)hold onezplot(cos(t),-pi/4,pi);grid on hold off例3：（1）（H:毕业设计程序niudun.m）function s=niudun(x,y,t)syms p;s=y(1);xishu=0;dxs=1;n=length(x);for(i=1:n-1) for(j=i+1:n) xishu(j) = (y(j)-y(i)/(x(j)-x(i); end temp1(i)=xishu(i+1); dxs=dxs*(p-x(i); s=s+temp1(i)*dxs; y=xishu;end simplify(s);if(nargin = 2) s=subs(s,p,x); s=collect(s); s=vpa(s,4); else m=length(t); for i=1:m temp(i)=subs(s,p,t(i); end s=temp;end(2)（H:毕业设计程序 niudun_main.m）x=1.0 1.3 1.6 1.9 2.2;y=0.7651977 0.6200860 0.4554022 0.2818186 0.1103623;yt=niudun(x,y,1.5)clear syms xfx=(63*x5-70*x3+15*x)/8;v=-1,1,-1,1;subplot(1,2,1);ezplot(fx),grid ontitle(5阶勒让德多项式)axis(v)x0=-1 -0.3 0.4 -0.7 0 -0.4 0.8 0.7 1;y0=subs(fx,x0);yt=niudun(x0,y0);subplot(1,2,2);ezplot(yt),axis(v),grid title(插值效果)t=0.24 -0.46 0.83legen5=subs(fx,t)chazhi=niudun(x0,y0,t)wucha=legen5-chazhi 例4：（H:毕业设计程序 Newtonforward.m与Newtonback.m）（1）function f= Newtonforward(x,y,x0)syms t; if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(x和y的维数不相等); return;end f = y(1);y1 = 0;xx =linspace(x(1),x(n),n);m1=(xx(2:n-1)= x(2:n-1);m=max(m1);if(m=1) disp(节点之间不是等距的); return;end for(i=1:n-1) for(j=1:n-i) y1(j) = y(j+1)-y(j); end c(i) = y1(1); l = t; for(k=1:i-1) l = l*(t-k); end; f = f + c(i)*l/factorial(i); simplify(f); y = y1; if(i=n-1) if(nargin = 3) f = subs(f,t,(x0-x(1)/(x(2)-x(1); else f = collect(f); f = vpa(f, 6); end endend（2）function f = Newtonback(x,y,x0)syms t; if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp(x和y的维数不相等); return;end f = y(n);y1 = 0;xx =linspace(x(1),x(n),n);m1=(xx(2:n-1)= x(2:n-1);m=max(m1);if(m=1) disp(节点之间不是等距的); return;end for(i=1:n-1) for(j=i+1:n) y1(j) = y(j)-y(j-1); end c(i) = y1(n); l = t; for(k=1:i-1) l = l*(t+k); end; f = f + c(i)*l/factorial(i); simplify(f); y = y1; if(i=n-1) if(nargin = 3) f = subs(f,t,(x0-x(n)/(x(2)-x(1); else f = collect(f); f = vpa(f, 6); end endEnd例5：（H:毕业设计程序herm.m与herm_main.m）(1)hermite插值的符合计算版%如果 输入向量只有三个 没有最后要插值的点%那么函数直接输出多项式%如果包含插值点，那么给出插值结果%其中要插值的点，可以是一个数，一个