1、第三章 插值法与最小二乘法 在实际问题中碰到函数 有些有解析表示式，但很复杂，有些只给出一些离散数据，这给我们求解函数值、导数值、零点、极值和积分值等带来了很多不便。对于这些情况，自然想法是，设法找到某个简单近似函数满足 。本章介绍两种方法，即插值法和最小二乘法。1/573.1 插值法 3.1.1 多项式插值概念在众多插值函数中，多项式是最简单最易计算。设 函数在区间a,b上连续，且在n+1个不一样点上分别取值为 。在多项式插值中，最基本、最简单问题是求一个次数不超出 代数多项式（3-1）满足其中 均为实数（3-2）2/57称为被插值函数；称插值多项式；条件（3-2）为插值条件；为插值点。其几
2、何意义为图3-1 插值几何意义满足（3-2）式插值多项式是存在且唯一。原因是满足3/57（3-2）式多项式未知系数行列式为著名范德蒙德（Vardermonde）行列式故解是存在且唯一。3.1.2 拉格朗日插值多项式 在每个插值点 结构插值基函数 为一个 次多项式 ，且满足条件 （3-3）4/57（3-4）由（3-4）式和多项式 定义以及（3-2）式插值条件，我们有：再由（3-4）式知，是 零点，故按因式定理必含有因式：既是 次多项式，可知 5/57再由条件 ，得于是有所以拉格朗日插值多项式可写为（3-5）6/57例3-1 已知特殊角 ，用 近似值。解：令用 为节点，有一次和二次拉格朗日插值公式
3、求7/57用 为节点，有8/57拉格朗日插值多项式程序设计计算公式(3-5)程序为二重循环。由内循环（j循环）经过累乘求得基函数 ，然后经过外循环（i循环）得到插值结果y.下列图描述了拉格朗日算法流程，对于该框图中相关参数说明以下：N:插值多项式次数（插值点个数-1）Xi：插值点Yi:插值点上函数值LN:拉格朗日插值结果T:存放累乘积I:外循环变量J:内循环变量9/57STARTINPUT N,X,Xi,Yi(I=0,N)LN=0I=0,NT=1J=0,NJ=I?NOT=(X-Xj)/(Xi-Xj)*TEND JLN=LN+Yi*TEND IOUTPUTENDYES10/57拉格朗日算法程序D
4、IMENSION X(I),Y(I)REAL LNREAD(*,*)NDO K=1,N READ(*,*)X(I),Y(I)END DOLN=0.0DO I=0,N T=1 DO J=0,N T=T*(X-X(J)/(X(I)-X(J)END DO LN=LN+T*Y(I)END DOWRITE(*,*)LN=,LN11/571、插值余项对插值后，用代替必定会产生误差，其误差可表示为 称上式为插值多项式 插值余项插值余项。设在区间上有直到阶导数，为 上个互异节点，为满足条件 次插值多项式。那么对于任何，有 (3-6)其中 且依赖于 。（普通取 最大值 3.1.3 插值多项式余项12/573.1
5、.4牛顿插值多项式 在插值基函数中，按这种方式来求：由线性代数性质可知，一个不高于n次多项式可表示成线性组合。在满足插值条件情况下，可表示为（记为）：（3-6）上式称为牛顿（牛顿（Newton）插值多项式）插值多项式。13/573.1.4牛顿插值多项式1.1、向前差分与牛顿向前插值多项式 取节点为等距离，则 其中h为步长。在两相临节点处函数值之差为 称为函数f(x)在节点处以h为步长一一阶阶向向前前差差分分（简称一阶差分）。又称一阶差分差分为二阶差分，记为 14/573.1.4牛顿插值多项式为了便于计算与一目了然，用表格描述为 表3.115/573.1.4牛顿插值多项式在等距离节点情况下，能够
6、用差分来确定(3.6)式待定系数，并将牛顿插值多项式加以简化。在节点上，函数值是已知，所以 16/573.1.4牛顿插值多项式将上式带入(3.6)式，并令 ，则式(3.6)可简化为 该式称为牛顿向前插值公式牛顿向前插值公式。其插值余项是（3-7）（3-8）17/573.1.4牛顿插值多项式2、向后差分与牛顿向后插值公式a、向后差分b、中心差分 18/573.1.4牛顿插值多项式3、差商与牛顿插值多项式用差商来确定(3.6)式中系数一阶差商 二阶差商 依次类推 19/573.1.4牛顿插值多项式表3.差商表20/573.1.4牛顿插值多项式x4.00024.01044.02334.0294f(x
7、)0.6020817 0.6031877 0.6045824 0.6052404例3-2 给出列表函数f(x)=lgx值，以下表所表示，试用牛顿插值法求lg4.01值。表3.3 解：依据给定函数作差商表，如表3.3所表示。x=4.01,x1=4.0002,x2=4.0233,x3=4.029421/573.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.604582
8、40.60524040.1084330.1081160.107869-0.013636-0.0130000.021781f(x)一阶差商二阶差商三阶差商将以上各插值点值、函数值和差商值代入插值公式得lg4.01=0.603144322/573.1.5 Hermite插值多项式 前述插值多项式均只要求插值点函数值相等，为深入光滑函数曲线，有时还要求其导数值也相等，并等于已知值 次多项式 。该多项式 称为 Hermite插值多项式，它也是一个近似式。求解Hermite插值多项式一个最简单方法就是直接应用牛顿插值法。即已知节点 处函数值与导数值 时，把 视为 重节点 ，并注意 个重节点 阶差商23/
9、573.1.5 Hermite插值多项式例 3-2 已知函数 函数值、导数值以下，求其插值多项式及误差表示式。-1010-4-2056解：24/573.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-1102125/573.1.6 三次样条插值多项式 1、三次样条插值函数定义 对于一组已知数据 ，且 ，若函数 满足：(1)、在每一个子区间上都能够用最高为三次多项式来表示。、在上二次连续可微。即，连续。、插值条件
10、 则称为函数f(x)关于节点 三次样条插值函数三次样条插值函数。26/573.1.6 三次样条插值多项式 2、边界条件问题 依据定义，在任意子区间上，三次样条函数可表示为 该式中有4个待定系数，故共有4n个待定系数。由定义中连续条件和插值条件可列出以下方程：27/573.1.6 三次样条插值多项式上式共有4n-2 个方程，所以还需要2个方程才能确定4n个系数，这就要用到边界条件。边界条件类型很多，常见有三种：1、在边界上给定一阶导数值；2、给定二阶导数值；（假如，称其为自然边界条件自然边界条件3、f(x)是以b-a为周期函数，则要求S(x)及其导数也是以b-a为周期函数，对应边界条件为 28/
11、573.1.6 三次样条插值多项式3、三弯矩方程设S(x)在某一节点二阶导数为（在力学上解释为：Mi是细梁在截面处弯矩，因为该弯矩与相临两个弯矩相关，故称三弯矩方程）因为S(x)最高次数不超出三次，所以它二阶导数是线性函数。令，于是29/573.1.6 三次样条插值多项式连续积分两次得其中 为积分常数。由以下插值条件确定 所以（3-9）30/573.1.6 三次样条插值多项式代入(3.9)式得 从上式可知，只要确定M值（M0,M1,Mn 共n+1个），就能够完全确定S(x)。在确定(3.9)式中待定系数时用了插值点条件，下面利用节点一阶导数连续条件（3-10）（3-11）31/573.1.6
12、三次样条插值多项式由(3.10)式得于是（3-12）（3-13）32/573.1.6 三次样条插值多项式由(3.11)式、(3.12)式和(3.13)式得若记则方程可简写为（3-14）33/573.1.6 三次样条插值多项式即 上方程组有n-1个方程，但有n+1个未知数。所以，要确定Mi还需用边界条件。1、第一个边界问题 边界条件则 （3-15）34/573.1.6 三次样条插值多项式整理得 其中（3-16）35/573.1.6 三次样条插值多项式2、第二种边界问题边界条件 次边界条件表明M0和Mn是已知，所以(3.15)式直接改写为（3-17）36/573.1.6 三次样条插值多项式2、第二
13、种边界问题边界条件a、b、并注意到 ，得记37/573.1.6 三次样条插值多项式则方程可简写成将上式和 与(3.15)式合在一起 满足第1或第2或第3种边界条件三次样条插值函数S(x)是存在且唯一。（3-18）38/573.1.6 三次样条插值多项式例例3-33-3：已知函数y=f(x)函数值以下 x -1.5 0 1 2 y 0.125 -1 1 9 求三次样条函数S(x)，使其满足边界条件 。解：这是第一个边界条件问题，(3-14)式和(3-16)式个参数为 39/573.1.6 三次样条插值多项式方程组为 求解后得 由(3-10)式得 40/573.2 曲线拟合 3.2.1 问题描述试
14、验目标是寻找一些物理量之间关系，如一组试验数据（xi,yi）,要寻找函数一个近似表示式，这就是一个曲线拟合问题。插值方法能够在一定程度上处理曲线拟合问题，但有显著不足。首先，试验数据有误差，假如进行插值，则所求得曲线必须严格经过全部试验数据点，从而使得该曲线保留了试验误差；其次，试验数据往往很多，用插值法求得表示式所以缺乏实用性。曲线拟合还有一个好坏标准问题，标准不一样决定着不一样拟合方法。惯用有最小二乘法，其基本思想是：使拟合曲线与试验点之差平方和最小。41/573.2 曲线拟合3.2.2 概述 试验数据表 试验编号 自变量 因变量 1 x1 y1 2 x2 y2 3 x3 y3 i xi
15、yi m xm ym42/573.2 曲线拟合 设x和y之间为线性关系，则有：经过试验数据 来确定（3-19）式中常数a和b，即建立x与y关系，此为一个一元线性拟合问题（线性回归）。问题：怎样用m组数据确定线性方程中a和b？见下列图。（3-19）43/573.2 曲线拟合 怎样确定哪一条直线是最好？a.最小二乘法：使回归残差平方和最小。b.残差：试验数据yi与回归方程计算f(xi)之间差，用qi 表示：按最小二乘法定义：（3-20）（3-21）44/573.2 曲线拟合2.算法 由（3-21）式可知，Q是a和b函数，即依据数学知识，要使Q最小，有：,将（3-21）式代入前两式：45/573.2
16、 曲线拟合 经过计算得到：经过计算得到：（3-22）（3-23）46/573.2 曲线拟合由（3-22）和（3-23）式可验证：3.方差分析 怎样衡量回归直线与试验数据之间吻合程度？或回归直线可信度是多少？-这个问题由方差分析来完成。衡量标准：（几个统计量）a.残差平方和 b.回归平方和 其中 47/573.2 曲线拟合c.剩下标准差 其中n=1 d.相关系数e.综合检验值 普通惯用是：相关系数R（R:0=R=1;R越靠近1越好）和综合检验值F（F为显著性检验，它值越大越好）。48/573.2 曲线拟合v多元线性最小二乘法：自变量个数大于一，如：当Q/bi=0时，得到一个线性方程组，用以求出
17、bi v非线性最小二乘法：将非线性方程回归转化为线性形式。如 转化为 即相当于v若非线性方程不能回归转化为线性形式时，可采取后面求根法中迭代法。49/573.2 曲线拟合3.2.3程序设计程序功效是：由已知m组实测数据（xi,yi),i=1,2,3,m,按最小二程序原理拟合多项式系数。其中用户能够依据详细问题预先确定拟合多项式次数，也能够按精度要求，使程序自动选择多项式次数。应该注意是，不论何种选择，其次数都不应该超出m。1、变量及数组使用说明M:实测数据（xi,yi)个数；N:拟合多项式次数加1，例初值为2，即为一次多项式；EPS:允许误差；N0:用户直接选择多项式次数，=0时表示按精度要求
18、自动选择多项式次数；50/573.2 曲线拟合X(M):存放给定点值；Y(M):存放给定点实测值；F(M):存放给定点所拟合函数值；S(2M):存放运算中形成正规方程组系数全部值；A(M,M):存放正规方程组系数增广矩阵；B(M):存放正规方程组常系数列向量；Z(M):调用解方程组子程序得到解向量；C(M):存放多项式系数列向量；Q(M):存放残量之绝对值。2、程序流程图51/57开始输入 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调用主元法解方程
19、组子程序，解向量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+1NQ0)THEN510 Q0=Q(I)520 I0=I54/573.2 曲线拟合520 ENDIF 670 RETURN 530 END DO540 IF(N0=0)THEN550 IF(Q0/F(I0)=EPS)THEN570 N=N+1580 IF(N0,a,b为待定参数)，使它能和以下观察值拟合。x 1 2 3 4 5 6 7 8 y 15.3 20.5 27.4 36.6 49.1 65.6 87.8 117.62、已知函数有y=f(x)函数值以下 x 0 1 2 3 4 y -8 -7 0 19 56求三次样条插值函数S(x)，使满足边界条件：S(0)=0,S(4)=4856/57上机实习题画出牛顿插值多项式程序结构流程图并用一个高级语言编制对应程序计算函数f(x)近似牛顿插值多项式N5(x),x=0.2+0.01k,k=020 x 0.2 0.24 0.28 0.32 0.36 0.40 y 0.1987 0.2377 0.2764 0.3146 0.3523 0.3894 57/57
