插值法与最小二乘法课件

插值法与最小二乘法

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

§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.5

Hermite插值多项式

前述的插值多项式均只要求插值点函数值相等，为进一步光滑函数曲线，有时还要求其导数值也相等，并等于已知值的次多项式。该多项式称为的Hermite插值多项式，它也是一种近似式。求解Hermite插值多项式的一种最简单方法就是直接应用牛顿插值法。即已知节点处函数值与导数值时，把视为重节点，并注意个重节点的阶差商3.1.5

Hermite插值多项式例3-2已知函数

的函数值、导数值如下，求其插值多项式及误差表达式。-1010-4-2056解：3.1.5

Hermite插值多项式作差商表如下。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概述实验数据表

实验编号

自变量因变量

1x1

y12x2

y23x3

y3

…

…

…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、源程序清单10

DIMENSIONX(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)270ENDDO

280ENDDO290C

产生增广矩阵300DOI=1,N310B(I)=0.0320DOJ=1,M330B(I)=B(I)+Y(J)*X(J)**(I-1)340ENDDO350A(I,N+1)=B