数值分析课程报告

插值法和多项式拟合的研究摘要在科研和生产实践中，常常需要通过一组测量数据来寻找变量x与y的函数关系近似表达式。解决这类问题的方法有两种：一种是插值法，另一种是拟合法。插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。拟合法能够是从给定的一组实验数据出发，寻找函数的一个近似表达式，该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合。本文主要介绍拉格朗日插值法、埃尔米特插值法、三次样条插值法以及基于最小二乘法的多项式拟合。关键词：拉格朗日插值，埃尔米特插值，样条插值，多项式拟合1方法的意义在许多实际问题及科学研究中，因素之间往往存在着函数关系，然而，这种关系经常很难有明显的解析表达，通常只是由观察与测试得到一些离散数值。有时，即使给出了解析表达式，却由于表达式过于复杂，不仅使用不便，而且不易于进行计算与理论分析。解决这类问题的方法有两种：一种是插值法，另一种是拟合法。插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。它要求给出函数的一个函数表，然后选定一种简单的函数形式，比如多项式、分段线性函数及三角多项式等，通过已知的函数表来确定一个简单的函数(X)作为f(x)的近似，概括地说，就是用简单函数为离散数组建立连续模型。插值法在实际应用中非常广泛，但是它也有明显的缺陷，一是测量数据常常带有测试误差，而插值多项式又通过所有给出的点，这样就是插值多项式保留了这些误差；二是如果实际得到的数据过多，则必然得到次数较高的插值多项式，这样近似的效果并不理想。拟合法能够很好的解决这些问题，它从给定的一组实验数据出发，寻找函数的一个近似表达式y=「(x),该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合的问题，函数的近似表达式丫=“)称为拟合曲线。常用最小而二乘法来确定拟合曲线。2插值法的介绍2.1插值法定义设f(x)为［a,b］上的函数，在互异点xo,Xp…,xn处的函数值分别为f(x。)，f(xj,・・・，f(Xn),构造一个简单函数(x)作为函数f(x)的近似表达式y=f(x)」(x)，使g二f(X),i=0,1,2，…,n(1.0)则称(X)为关于节点X0,X-!,・・・,Xn的插值函数；称x0,x1,・・・,xn为插值节点；称(Xif(Xi)),i=1,2，…，n为插值点；f(x)称为被插值函数。式(1・0)称为插值条件。这类问题称为插值问题。插值的任务就是由已知的观测点，为物理量(未知量)建立一个简单的、连续的解析模型，以便能根据该模型推测该物理量在非观测点处的特性。常用的插值函数类｛「(X)｝是代数多项式，相应插值问题是代数插值，本文主要介绍三种代数差值：拉格朗日插值，埃尔米特插值和样条插值。2.2拉格朗日插值2.2.1两点插值问题已知y=f(x)(=0,1)，求满足插值条件R(x)=y(i=o,i)的插值多项式R(x)。根据解析几何知识可知，所求的R(x)为过点(xo,yo),(Xj,y)的直线，即:R(x)=yXo+—(x-xo)i—xox-Xo显然lo(Xi)=(i=O)X.-Xiooio(1=1)1(i=1)上式经整理可改写为：式中l°(X)二li(x)且lo(X)，h(X)为由插值节点唯一确定的线性函数。l°(x),l’(x)为节点Xo,Xir上的一次插值基函数。可以看出，节点XO,X1上的插值基函数的次数为插值节点个数减一，基函数组中所含的函数个数与插值节点数相同。而满足R(Xi)=yi(i=O,1)的插值多项式P|(x)就是节点xo，x1上插值基函数的线性组合，其组合系数分别为yo,y1。这种表示为插值基函数线性组合的一次插值多项式也就是一次拉格朗日插值多项式。当给定n+1个插值节点后，可类似定义n次插值基函数，并以此构造n次拉格朗日多项式。2.2.2n次拉格朗日插值多项式n次插值基函数：lk(x)二(k=O,1...n)(X_X。)…(X—Xkj(Xlk(x)二(k=O,1...n)(Xk—xo)・・・(XkXkJ(Xk-Xk1)…(Xk-Xn)显然lk(x)具有以下性质:性质1，lk(x)」_(k=O,1...,n)O(i=k)性质2,lk(x)(k=O,1...,n)为由插值节点Xo,X15・・・,Xn唯一确定的n次函数性质3,基函数组所含的基函数个数与插值节点个数相同

以n次插值基函数为基础，可得拉格朗日插值多项式为：L(X)-l/X)\,L(X1Xo)(X1Xk」)(X1Xk"...(X1Xn)yL(X)lX)o)…kknyLn(X)lk(X)yknk£7%—X0)...(Xk—x/Xk—兀出)…区一Xn(nn八(i【上乞皿k卫」兀-乂TOC\o"1-5"\h\zkiji#nn若记「nl(X)(X-Xj),则有51(XQ八(兀-人)。于是-(X)也可写成:7ii=#0nl(X)Ln(X)nl心(X—Xk>3n+(Xk)223拉格朗日插值多项式的余项Rn(X)=f(X)-Ln(X)称为拉格朗日插值多项式的余项，也叫截断误差。yk用简单的插值多项式Ln(x)代替复杂的函数f(x)yk用简差是否满足精度要求。拉格朗日型余项：RM二f(X)-Ln(X)(]((X)(n+1)!Rn(x)(n1)33(x)其中宀1&)「(x-X)，'-(a,b)且依赖于X。。—般情况下,匚e(a,b(的具体i=Q数值无法知道，但是若能够求出maxf(ni)(x)Rn(x)(n1)33(x)由此看出，Rn(x)的大小除了与Mm有关外，还与插值节点有密切关系。当给定m个点出的函数值，但仅选用其中n1(nT:::m)个作为插值条件而求某点x处函数值时，n+1个节点x0,x15・・・,xn的选取应该尽可能的接近x,使得计算的函数值的误差限尽可能的小。2.3埃尔米特插值

许多实际问题不但要求插值多项式与被插值函数在节点处的函数值相当，而且还需要求其导数值相等。满足这种要求的插值多项式就是埃尔米特插值多项式。—般情形的埃尔米特插值问题—般情形的埃尔米特插值问题是指所满足的插值条件中函数值的个数与导数值的个数相等。即当函数f(x)在区间［a,b］上n+1个节点Xj(i=0,1，…,n)处的函数值yj二f(xi)及导数值mAf(xi)给定时，要求一个次数不超过2n+1的多项式H2ni(x),使之满足:日勿叩=久(i=0,1,・・・，n)H2nI")二m这里给出个2n+2个插值条件，可唯一确定一个形式为也川(刈二a。•ajX•a^；X2」1的多项式，但是要确定2n+2个系数非常复杂,因此此时可以借用构造拉格朗日插值多项式的基函数。设:“(xh/x)(j-0,1,...,n)为次数不超过2n+1的多项式，且满足：(i,j=0,1,・・・,n)：.(xj='••,匸.(X.)=0-(i,j=0,1,・・・,n)•(XAO,■(G•■J-j(G-j则满足上述条件的埃尔米特插值多项式可以写成用插值基函数表示的形式:nH2n1j(x)y.：j(X)mj]上式满足插值条件。可以得到基函数:j(X),'j(X)的解析式为:_1—_1—2(x_X/为Ik=0Xj-兀lj2(x)(J=0,1,…n),Pj(x)=(x—Xj)lj2(x)(j=0,1,..n),式中lj(x)为拉格朗日插值基函数因此埃尔米特插值多项式即为：n|n|H2n1凶八1_2(x_Xj)Ej=0lj2(x)y.、(冷勺^网■j=般情形下的埃尔米特插值多项式的余项为:(2n+2)!式中：；：=(a,b),且与x有关。埃尔米特插值的几何意义：曲线y二tni(x)与曲线y=f(x)在插值节点处有相同的公共切线。在带导数的插值问题中,有时插值条件中的函数值个数与导数值个数不等。这时可以以一般情况的埃尔米特插值多项式为基础,运用待定系数法求出满足插值条件的多项式。2.4样条插值在上述方法中，我们根据区间［a,b］上给出的节点得以得到函数f(x)的插值多项式，但是并非插值多项式的次数越高，逼近函数f(x)的精度越好，主要原因是因为对于任意的插值节点，当n》时，插值多项式Pn(x)不一定收敛到f(x)。这种高次插值不准确的现象称为龙格现象。为了避免高次插值的缺点，人们常常采用分段插值的方法,即将插值区间分为若干个小区间,在每个小区间上运用前面介绍的插值方法构造低次插值多项式。采用分段现象插值与分段二次插值,可以构造一个整个连续的函数,而采用分段三次埃尔米特插值则可以构造一个整体上具有一节连续导数的插值函数。实际问题中,很少给出插值点上的导数值。三次样条插值就是在只给出插值点上函数值的情况下,构造一个整体上具有二阶连续导数的插值函数。2.4.1三次样条插值函数的定义设区间［a,b］上取n+1个节点a=X。；咅：：：x2::…：：：xnJ<x,d=b，若函数S(x)满足条件：在整个区间［a,b］上具有二阶连续导数⑵在每个小区间［Xy，x］(i=0,1，…,n)上是x的三次多项式(3)S(x)=yi(i=0,1,••n),则称S(x)为f(x)的三次样条插值。S(x)的边界条件为：给定两端点处的一阶导数值,记为：S(Xo)=m°，S(Xn)=m给定两端点处的二阶导数值，记为:S(X°)=M°,s(Xn)二M此外，对S(Xo)=S(Xn)=0的边界条件称为自然边界条件2.4.2三次样条插值函数的构造构造三次样条插值函数，就是要写出它在子区间［X」xj(i=0,1n)上的表达式，记为S(x)(i=0,1，…,n)。一、用节点处一阶导数表示的三次样条插值函数记节点处的一阶导数值为S(x)二mi(i=0,1,...,n)，若已知mi后，则S(x)在［X」X］(i=0,1，…,n)上就是满足条件S(xi」)=yi」，S(xi)=yi，(i=0,1,...,n)S(4)=m」，S(x)=mi(i=0,1,-5n)的三次埃尔米特插值多项式。步骤如下：根据S(x)，S(x)在内阶段的连续性及插值条件，运用［Xj/Xj上的二点三次埃尔米特插值多项式，写出S(x)用mi(i=0,1,...5n)表示的形式。利用S(x)在内节点Xi(i=0,1，…，n-1)的连续性及边界条件，导出含mi(i=0,1,...,n)的n+1阶线性方程组。求解含mi(i=0,1,...,n)的线性方程组，将得到的m：代入伙二必］的二点三次埃尔米特插值多项式具体公式为：Si(x)=1+22Si(x)=1+22X—XiY,X八XXi-Xii—Xiji4XXi4i4+仆2x—IX—XjLXi^-Xi+(X-Xi)X一Xi!mu+(x—xjX—诊mi1g—X丿Xi—Xi斗丿1mi的确定由S(x)得二阶导数在内节点连续的条件得到，即：im|im|42miMm1二(i=0,1,..n,1)该式即为n+1阶线性方程组。其中：hi=Xi-X」hh1M1.i=1-■j=3（Mjh（八hi=Xi-X」hh1M1.i=1-■j=3（Mjh（八0,1,•叩）'丿经)hi申hi4if21扎2八2moM12

++2m22M4n4」12■2M1*-22M2+++》-n24n+1对于第二种边界条件，将其代入，得到其中：n_22m1>m2mA阶线性方程组:/rr,f2*-Mn/mnJn-n-1mn求解该线性方程组还需两个边界条件。对于第一种边界条件，将其代入，得到只含n-1个未知数的线性方程组:fo=3坐必-九0h12fn=3上归hn二、用节点处二阶导数表示的三次样条插值函数记节点处的二阶导数值为S(xJ=m,i=0,1,...,n)，由于S(x)在［八4旳(i=0,1，…,n)上是x的线性函数，因此构造以节点处二阶导数表示的三次埃尔米特插值多项式步骤如下：根据S(x)在内节点的连续性及为线性函数的特点，将S(x)表示为线性函数。在根据S(x)在内节点的连续性及插值条件，写出S(x)用M』=0,1,...,n)表示的形式。利用S(x)在内节点xi(i=0,1,...,n-1)的连续性及边界条件，导出XX—II求解含M/i=0,1,...,n)的线性方程组，将得到的Mi代入［Xj,X」的S(x)的表达式，即得到二点三次埃尔米特插值多项式具体公式为：S.(x)=AAMi」AXhAMi(%」Mj」十、为一xMiiiiii6h6hi-rh)〒(yriMi的确定由S(x)得二阶导数在内节点连续的条件得到，即：MM.2M「：；'％Mi4£(i=0,1,..n,1)iiji4-£该式即为n+1阶线性方程组。其中:hi=xh1hi1(i71,..n厂hi1(i71,..n厂1)i6hihidii1(片〔i1-yihi1解该方程还需要两个边界条件。对于第一种边界条件，将其代入，得到只含解该方程还需要两个边界条件。对于第一种边界条件，将其代入，得到只含_21【M12人+++Mn42打4式中：n-1个未知数的线性方程组「m1irf0'TOC\o"1-5"\h\z(iiiM2Ifl*I■■一1I1rpiM心|fz丿Ifn」X-■xy.jf对于第二种边界条件，将其代入,6(Xoh|h|-mo)6/(mn-hn得到Vn-yn4)hnn+1阶线性方程组:，，zM1'm2■FfM1M0f2aM22n-2M'n-22M2n2AM,nJ对于三次样条插值函数来说，n4.nj当插值节点逐渐加密时,fn.(fn丄—扎n/M可以证明：不但样条插值函数收敛于函数本身，而且其导数也收敛于函数的导数。3.曲线拟合3.1用最小二乘法求解矛盾方程组工程实际中许多问题都可归结为矛盾方程组，实际中需要寻求矛盾方程组的N未知数Xj,X2,・・・,xn的一组取值，它使得偏差的绝对值之和送冋尽可能的小，为了7NNn2便于分析计算和应用，常采用使偏差的平方和Q=vPjXjbJ达到i=4j=4最小值，这一条件称为最小二乘原则。按最小二乘原则来选择X4,X2,・・・,Xn的一组取值的方法称为求解矛盾方程组的最小二乘法，符合条件的人，X2,・・・,Xn的一组取值称为矛盾方程组的最小二乘解。把Q看出n个自变量Xj,X2,・・・,Xn的二次函数，记为Q=f(XpX2,・・・,Xn),因此，n求矛盾方程组的最小二乘解就是二乘函数的最小值点。若矛盾方程组-BjXAbij吕的系数矩阵A的秩为n,则二次函数Q=f(Xi，X2，...,Xn)—定存在最小值。即矛盾方程组的最小二乘解存在，且正则方程组AtAx二ATb有唯一解，此解就是矛盾方程组的最小二乘解。3.2多项式拟合设通过测量得到函数y二f(x)的一组数据为：y二f(X)(i=0,1，…，n),求一个次数低于N-1的多项式y=(x)=a0•a4X-a2X2...■amXm(m：：：N-1)(其中a0,a4,・・・,am待定)，使其最好的拟合这组数据，最好的标准时：使得(x)在Xi的偏差「二「(X)r偏差「二「(X)r的平方和Q八f「NN(N-yj］达到最小。2a0,a4,・・・,am的解i#i=4i#i即正则方程组AtAx二ATb的解。NNNNa°N+ai送Xi+a2瓦x：+...+amZ小i-1i」\-1\-1NNNNN』a°送人七〔送x2+aAx;+...+am瓦xj*二送xj%i£i)i」i=i—aNNa。瓦录七正NXm++aAx亡+…m+amEx：m匹Xji£777-y由该正则方程组求得的唯一解代入拟合多项式y=「（x），即为所求。5.数值试验表1是1971年到1990年我国总人口的统计数字，试根据1971年到1985年这15年人口统计数字用下面几种方法预测未来20年的人口数字，并用图示的方法比较1986年到1990年间预测人口数字与实际统计数字的差异，在你所使用的几种预测方法中找出一组较为合理的预测方法。（1）指数形式y二aebx；（2）拉格朗日插值、埃尔米特插值、样条插值；（3）三次多项式拟合；（4）四次多项式拟合表1人口统计数字年份1971197219731974197519761977197819791980人口8.52298.71778.92119.08599.24209.37179.49749.62599.75429.8705年份1981198219831984198519861987198819891990人口10.007210.165410.300810.435710.585110.750710.930011.102611.270411.4333数值试验结果:1■指数形式拟合图形为:2■三次多项式拟合的效果:19901938198619841982198019781976197419729951010.51111519708b——实际统计曲线19901938198619841982198019781976197419729951010.51111519708b*实际统计值三次多项式拟合

3.四次多项式拟合的效果199019881986实际统计曲线卄实际统计曲线卄实际统计值四初多项式拟合19301976197L41972—-1972—-参考文献王正林，龚纯，等•精通MATLAB科学计算•北京：电子工业出版社，2009:135—157任玉杰•数值分析及其MATLAB实现]M]•北京：高等教育出版社，2007：584—642李信真，车刚明，等•计算方法[M]•西安：西北工业大学出版社，2010:150—167金一庆，陈越•数值方法[M]•北京：机械工业出版社，2000:165—176附录1■按指数形式拟合的程序：x=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851,10.7507,10.9300,1;y=[1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,