Cha插值法与最小二乘法

第三章插值法与最小二乘法在实际问题中遇到的函数有些有解析表达式，但很复杂，有些只给出一些离散数据，这给我们求解函数值、导数值、零点、极值和积分值等带来了诸多不便。对于这些情况，自然的想法是，设法找到某个简单近似函数满足。本章介绍两种方法，即插值法和最小二乘法。

§3.1插值法3.1.1多项式的插值概念在众多的插值函数中，多项式是最简单最易计算的。设函数在区间[a,b]上连续，且在n+1个不同的点上分别取值为。在多项式插值中，最基本、最简单的问题是求一个次数不超过

的代数多项式

（3-1）满足其中均为实数（3-2）称为被插值函数；称插值多项式；条件（3-2）为插值条件；为插值点。其几何意义为图3-1插值几何意义满足（3-2）式的插值多项式是存在且唯一的。原因是满足（3-2）式的多项式的未知系数行列式为著名的范德蒙德（Vardermonde）行列式故解是存在且唯一的。3.1.2拉格朗日插值多项式

在每个插值点构造插值基函数为一个

次多项式，且满足条件

（3-3）（3-4）由（3-4）式和多项式的定义以及（3-2）式的插值条件，我们有：再由（3-4）式知，是的零点，故按因式定理必含有因式：既是次多项式，可知

再由条件，得于是有因此拉格朗日插值多项式可写为

（3-5）例3-1已知特殊角，用

的近似值。解：令用为节点，有一次和二次拉格朗日插值公式求用为节点，有拉格朗日插值多项式程序设计计算公式(3-5)的程序为二重循环。由内循环（j循环）通过累乘求得基函数，然后通过外循环（i循环）得到插值结果y.下图描述了拉格朗日算法流程，对于该框图中的有关参数说明如下：N:插值多项式次数（插值点个数-1）Xi：插值点Yi:插值点上的函数值LN:拉格朗日插值结果T:存放累乘积I:外循环变量J:内循环变量STARTINPUTN,X,Xi,Yi(I=0,…N)LN=0I=0,NT=1J=0,NJ=I?NOT=(X-Xj)/(Xi-Xj)*TENDJLN=LN+Yi*TENDIOUTPUTENDYES拉格朗日算算法程序DIMENSIONX(I),Y(I)REALLNREAD(*,*)NDOK=1,NREAD(*,*)X(I),Y(I)ENDDOLN=0.0DOI=0,NT=1DOJ=0,NT=T*(X-X(J))/(X(I)-X(J))ENDDOLN=LN+T*Y(I)ENDDOWRITE(*,*)‘LN=‘,LN1、插值余项对插值后，，用代代替必必定会产生生误差，其其误差可表表示为称上式为插插值多项式式的插值余项。设在在区间间上上有直到到阶阶导导数，为为上个个互异异的节点，，为为满足条件件的次插值多项项式。那么对于任任何，有(3-6)其中且且依赖赖于。。（（一般取最大值3.1.3插值多项式式余项3.1.4牛顿插值多多项式在插值基函函数中中，按按这种方式式来求：：由线线性代数性性质可知，，一个不高高于n次的的多项式可可表示成的线性组合合。在满足足插值条件件情情况下，可可表表示为（记记为））：（3-6））上式称为牛顿（Newton）插值多项项式。3.1.4牛顿插值多多项式1.1、、向前差分与与牛顿向前前插值多项项式取节点为等等距离，则则其中h为步长。在在两相临节节点处的函函数值之差差为称为函数f(x)在节点处以以h为步长的一阶向前差差分（简称一阶阶差分）。。又称一阶阶差分的差差分为二阶阶差分，记记为3.1.4牛顿插值多多项式为了便于计计算与一目目了然，用用表格描述述为表3.13.1.4牛顿插值多多项式在等距离节点点的情况下下，可以用用差分来确确定(3.6)式的待定系系数，并将将牛顿插值值多项式加以以简化。在节点上，，函数值是是已知的，，所以3.1.4牛顿插值多多项式将上式带入入(3.6)式，并令，，则式(3.6)可简化为该式称为牛顿向前插插值公式。其插值余余项是（3-7））（3-8））3.1.4牛顿插值多多项式2、向后差分与与牛顿向后后插值公式式a、向后差分分b、中心差分分3.1.4牛顿插值多多项式3、差商与牛顿顿插值多项项式用差商来确确定(3.6)式中的系数一阶差商二阶差商依次类推3.1.4牛顿插值多多项式表3.２差差商表3.1.4牛顿插值多多项式x4.00024.01044.02334.0294f(x)0.60208170.60318770.60458240.6052404例3-2给出列表函函数f(x)=lgx的值，如下下表所示，，试用牛顿顿插值法求求lg4.01的值。表3.3解：根据给定的的函数作差差商表，如如表3.3所示。x=4.01,x1=4.0002,x2=4.0233,x3=4.02943.1.4牛顿插值多多项式lgx=f(x0)+f(x0,x1)(x-x0)+f(x0,x1,x2)(x-x0)(x-x1)+f(x0,x1,x2,x3)(x-x0)(x-x1)(x-x2)表3.4x4.00024.01044.02334.02940.60208170.60318770.60458240.60524040.1084330.1081160.107869-0.013636-0.0130000.021781f(x)一阶差商二阶差商三阶差商将以上各插插值点值、、函数值和和差商值代代入插值公公式得lg4.01=0.60314433.1.5Hermite插值值多项式前述的插值值多项式均均只要求插插值点函数数值相等，，为进一步步光滑函数数曲线，有有时还要求求其导数值值也相等，，并等于已已知值的次次多多项式。。该该多项式称称为的的Hermite插值多多项式，它它也是一种种近似式。。求解Hermite插值多项项式的一种种最简单方方法就是直直接应用牛牛顿插值法法。即已知知节点处处函函数值与导导数值时，把视视为重重节点，并注意个个重节点点的阶阶差差商3.1.5Hermite插值值多项式例3-2已知函数的函数值、、导数值如如下，求其其插值多项项式及误差差表达式。。-1010-4-2056解：3.1.5Hermite插值值多项式作差商表如如下。f(0,0)=y’(0)/1!=0,f(0,0,0)=y””(0)/2!=3,f(1,1)=y’(1)/1!=5.按差商商表3.2有：-1000110-4-4-4-2-2-400254323-1-110213.1.6三次样条插插值多项式式1、三次样条插插值函数的的定义对于一组已已知的数据据，，且且，，若函数满满足：：(1)、在每一个子子区间上上都都可以用最最高为三次的的多项式来表表示。⑵、在在上上的二次次连续可微。。即，，，，连续。⑶、插值条件件则称为函数f(x)关于于节点的的三次样条插值值函数。3.1.6三次样条插值值多项式2、边界条件问题题根据定义，在在任意子区间间上，三次样样条函数可表表示为该式中有4个待定系数，，故共有4n个待定系数。。由定义中的的连续条件和和插值条件可可列出下列方程程：3.1.6三次样条插值值多项式上式共有4n-2个方程，因此此还需要2个方程才能确确定4n个系数，这就就要用到边界界条件。边界条条件的类型很很多，常见的的有三种：[1]、在边边界上给定一一阶导数的值值；；[2]、给定定二阶导数的的值；（（如果，称其为自然边界条件件[3]、f(x)是以以b-a为周周期的函数，，则要求S(x)及其导导数也是以b-a为周期期的函数，相相应边界条件件为3.1.6三次样条插值值多项式3、三弯矩方方程设S(x)在在某一节点的的二阶导数为为（在力学上解解释为：Mi是细梁在截面面处的弯矩，，因为该弯矩矩与相临的两两个弯矩有关关，故称三弯弯矩方程）因为S(x)最高次数不不超过三次，，所以它的二二阶导数是线线性函数。令，，于是是3.1.6三次样条插值值多项式连续积分两次次得其中为为积分常数。。由以下插值值条件确定所以（3-9）3.1.6三次样条插值值多项式代入(3.9)式得从上式可知，，只要确定M值（M0,M1,…,Mn共n+1个）），就可以完完全确定S(x)。在确确定(3.9)式中待定定系数时用了了插值点的条条件，下面利利用节点一阶阶导数连续条条件（3-10））（3-11））3.1.6三次样条插值值多项式由(3.10)式得于是（3-12））（3-13））3.1.6三次样条插值值多项式由(3.11)式、(3.12)式和(3.13)式得得若记则方程可简写写为（3-14））3.1.6三次样条插值值多项式即上方程组有n-1个方程程，但有n+1个未知数数。因此，要要确定Mi还需用边界条条件。[1]、第一种边界问问题边界条件则（3-15））3.1.6三次样条插值值多项式整理得其中（3-16））3.1.6三次样条插值值多项式[2]、第二二种边界问题题边界条件次边界条件表表明M0和Mn是已知的，所所以(3.15)式直接改写为为（3-17））3.1.6三次样条插值值多项式[2]、第二二种边界问题题边界条件a、b、并注意到，，得记3.1.6三次样条插值值多项式则方程可简写写成将上式和与与(3.15)式合在一起满足第[1]或第[2]或第[3]种边界条件件的三次样条条插值函数S(x)是存在且唯一一的。（3-18））3.1.6三次样条插值值多项式例3-3：已知函数y=f(x)的函数值如下下x-1.5012y0.125-119求三次样条函函数S(x)，使其满足边边界条件。。解：这是第一一种边界条件件问题，(3-14)式式和(3-16)式的个个参数为3.1.6三次样条插值值多项式方程组为求解后得由(3-10)式得3.2曲线拟合3.2.1问题描述实验的目的是是寻找一些物物理量之间的的关系，如一一组实验数据据（xi,yi）,要寻找函数的的一个近近似表达式，，这就是是一个曲线拟拟合问题。插值方法可以以在一定程度度上解决曲线线拟合问题，，但有明显的的不足。首先先，实验数据据有误差，如如果进行插值值，则所求得得的曲线必须须严格通过所所有实验数据据点，从而使使得该曲线保保留了实验误误差；其次，，实验数据往往往很多，用用插值法求得得的表达式因因此缺乏实用用性。曲线拟合还有有一个好坏的的标准问题，，标准不同决决定着不同拟拟合的方法。。常用的有最最小二乘法，，其基本思想想是：使拟合合曲线与实验验点之差的平平方和最小。。3.2曲线拟合3.2.2概述实验数据表实验编号自变量因因变量量1x1y12x2y23x3y3……………ixiyi……………mxmym3.2曲线拟合设x和y之间间为线性关系系，则有：通过实验数据据来来确定（3-19）式式中的常数a和b，即建建立x与y的的关系，此为为一个一元线线性拟合问题题（线性回归归）。问题：如何用用m组数据确确定线性方程程中a和b？？见下图。（3-19））3.2曲线拟合如何确定哪一一条直线是最最好的？a.最小二乘法：：使回归的残残差平方和最最小。b.残差：实实验数据yi与回归方程计计算的f(xi)之间的差，，用qi表示：按最小二乘法法的定义：（3-20））（3-21））3.2曲线拟合2.`算法由（3-21）式可知，Q是a和b的函数，即根据数学知识识，要使Q最小，有：,,将（3-21）式代入前两两式：3.2曲线拟合经过计算得到到：（3-22））（3-23））3.2曲线拟合由（3-22）和（3-23）式可可验证：3.方差分分析如何衡量回归归直线与实验验数据之间的的吻合程度？？或回归直线线的可信度是是多少？----这个个问题由方差分析来完完成。衡量标准：（（几个统计量量）a.残差差平方和b.回归归平方和其其中3.2曲线拟合c.剩余标准差其其中中n=1d.相关系数数e.综合检验验值一般常用的是是：相关系数数R（R:0<=R<=1;R越接近1越好）和综合检验验值F（F为显著性检检验，它的值值越大越好）。3.2曲线拟合多元线性最小小二乘法：自自变量个数大大于一，如：：当Q/bi=0时，得到一个个线性方程组组，用以求出出bi非线性最小二二乘法：将非非线性方程回回归转化为线线性形式。如转转化为即即相当当于若非线性方程程不能回归转转化为线性形形式时，可采采用后面求根根法中的迭代代法。3.2曲线拟合3.2.3程程序设计程序功能是：：由已知的m组实测数据据（xi,yi),i=1,2,3,…,m,按最小二二程序原理拟拟合多项式的的系数。其中中用户可以根根据具体问题题预先确定拟拟合多项式的的次数，也可可以按精度要要求，使程序序自动选择多多项式的次数数。应该注意意的是，无论论何种选择，，其次数都不不应该超过m。1、变量及数数组使用说明明M:实实测数据据（xi,yi)的个数；N:拟拟合多项项式次数加1，例初值为为2，即为一一次多项式；；EPS:允允许许误差；N0:用用户直直接选择多项项式次数，=0时表示按按精度要求自自动选择多项项式次数；3.2曲线拟合X(M):存存放给定点的的值；Y(M):存存放给定点的的实测值；F(M):存存放给定点所所拟合的函数数值；S(2M):存放放运算中形成成的正规方程程组系数的所所有值；A(M,M):存放放正规方程组组系数增广矩矩阵；B(M):存存放正规方程程组常系数列列向量；Z(M):调调用解解方程组子程程序得到的解解向量；C(M):存存放放多项式系数数列向量；Q(M):存存放放残量之绝对对值。2、程序流程程图开始输入N0,M,X(I),Y(I)N0=0?N=2

N=N0+1计算S(K),K=1,2,…,2N-1生成S(I+J-1)=A(I,J)I=1,…N;J=1,…,N计算B(K)=A(K,N+1),K=1,…,N调用主元法解解方程组子程程序，解向量量Z(I)X(I)=Z(I)I=1,…,N计算拟合曲线线在给定点X(I)上函函数值F(I)I=1,…,M计算残量Q(I)=ABS(F(I)-Y(I))寻找最大值Q(I0)N0=0?Q(I0)/F(I0)<EPS?

N=N+1N<M输出多项式系系数C(I),I=1,…,M结束YNYNYNYN3.2曲线拟合3、源程序清清单10DIMENSIONX(M),Y(M),20&S(2*M),A(M,M),B(M),30&Z(M),P(M),Q(M)40C输入已知数据据50READN0,M60IF(N0=0)THEN70READEPS80ENDIF90DOI=1,M100READ(*,*)X(I),Y(I)110ENDDO120IF(N0=0)THEN130N=2140ELSE150N=N0+1160ENDIF170C生生成系数矩矩阵S(2N-1)180DOI=1,2*N-1190S(I)=0.0200DOJ=1,M210S(I)=S(I)+X(J)**(I-1)220ENDDO230ENDDO240DOI=1,N3.2曲线拟合250DOJ=1,N260A(I,J)=S(I+J-1)270ENDDO280ENDDO290C产生增广矩阵阵300DOI=1,N310B(I