《应用数值分析》课件数值分析2.3牛顿插值法.ppt

1、第1.3节 牛顿插值法差商Newton插值由由 组不同数据组不同数据 构造的构造的 次多项式次多项式1n n0(,)1,，iixy in L0011()()()()nnnLxlx ylx ylx yL12111211()()()()()()()()()()()kknkkkkkkkknxxxxxxxxxxlxxxxxxxxxxx LLLL其中其中拉格朗日插值多项式拉格朗日插值多项式0()nkkklx y 含义直观含义直观形式对称形式对称优点：优点：缺点：缺点：增加一个节点时，全部基函数增加一个节点时，全部基函数 都都需重新算过。需重新算过。公式不具有继承性，不利于编程。公式不具有继承性，不利于编

2、程。()klx0 01 1()()()()nn nLxy lxy lxy lxL分析：分析：在n次多项式空间中另找一组合适的基函数01(),(),()nxxxL0011()()()()nnnP xcxcxcxL希望：希望：每增加一个节点时，只需重算一个基函数Newton Newton 插值：插值：010201011()()()()()()()nnnN xcc x xc x xx xc x xx xx x LL系数 是多少？疑问：kc容易求吗？容易求吗？差商差商(亦称均差亦称均差)/*divided difference*/01201(),1,(),(,),(nnf xa bnx xxxf xf

3、 xf x函数在区间上个互异节点处的函数已值知:()(),(),jiijjiijf xf xf x xxf xxxx称为关于节：点的一阶差商,jkijijkkif xxf x xf x xxxx(),ijkf xx xx为关于节点的二阶差商110110110,kkkkkkf xxxf x xxf x xxxxx01(),kf xxxxk为关于节点的阶差商差 商 表f(x3)x3f(x2)x2f(x1)x1f(x0)x0f x1,x2,x3f x0,x1,x2二阶差商f x2,x3f x1,x2f x0,x1一阶差商f(x)xf x0,x1,x2,x3三阶差商由差商定义可知：高阶差商是两个低一阶

4、差商的差商。,jkijijkkif xxf x xf x xxxx()()()()kjjikjjikif xf xf xf xxxxxxx()()()()()()()()()kjjijikjkjjikif xf xxxf xf xxxxxxxxx()()()()()()()()()()()kjijjijkjikjkjjikif xxxf xxxf xxxf xxxxxxxxx()()()()()()()()()jkikjkikjjijikif xf xf xxxxxxxxxxxxx()()()()()()()()()jkikjkijkjiijikf xf xf xxxxxxxxxxxxx()(

5、)(),()()()()()()jikijkijikjkjikjkif xf xf xf x x xxxxxxxxxxxxx01()(),(),()kf xkf xf xf x的差商可以由函值性表示阶数线为(差商可用函数值的线性组合表示差商可用函数值的线性组合表示)：性质性质10110011(),()()()()kikkiiiiiiikf xf x xxxxxxxxxxx()()(),()()()()()()jikijkijikjkjikjkif xf xf xf x x xxxxxxxxxxxxx,ijx x差商中任意交的序，差商值不。换两个节点顺变(差商具有对称性差商具有对称性)：性质性质

6、2110201,=kkkkkf xxf xxxxxx差商值与节点顺序无关差商值与节点顺序无关01(),kf xa bkx xxa b若 在 上存在，且阶导数节点(差商与导数的关系差商与导数的关系)：性质性质3()011(),=!kkkff x xxxk,a b 则，有Newton Newton 插值：插值：010201011()()()()()()()nnnN xcc x xc x xx xc x xx xx x LL系数系数 是多少？是多少？疑问：疑问：kc容易求吗？容易求吗？Newton Newton 插值：插值：010201011()()()()()()()nnnN xcc x xc x

7、 xx xc x xx xx x LL (n)(n-1)(2)(1)(1)000()(),()f xf xf x x x x(2)001011,()f x xf x xf x x x x x010120122,()f x x xf x x xf x x x x x x(3)(n)0110101,()nnnnf x xxxf x xxf x xx x x x0010012010120110101()(),(),()(),()()(),()()()nnnnf xf xf x xxxf x x xxxxxf x x xxxxxxxxf x xx x xxxxxx0010012010120110101

8、()(),(),()(),()()(),()()()nnnnf xf xf x xxxf x x xxxxxf x x xxxxxxxxf x xx x xxxxxx是是n 次多项式次多项式满足插值条件满足插值条件001001201012011()(),(),()(),()()()nnnN xf xf x xxxf x x xxxxxf x x xxxxxxxx011(),()nnnR xf x xx xx(1)1()()()()(1)!nnnnfN xL xxn()011()3,=!kkkff x xxxk性质Newton Newton 插值：插值：010201011()()()()()()

9、()nnnN xcc x xc x xx xc x xx xx x LL011,(0,1,2,)kkkcf x xxxkn其 中，优点：优点：每增加一个节点，只要再增加一项即可每增加一个节点，只要再增加一项即可 即有递推公式即有递推公式101101()(),()()()nnnnnNxN xf x xx xxxxxxx例:已知函数已知函数 的数值表的数值表 sinx00030ix0045iy0.500000.707110600900.866031试求试求 的四次的四次Newton插值多项式，并计算插值多项式，并计算 的近似值。的近似值。sinx0sin40解：解：取节点取节点 ,则则 的差商表如

10、下的差商表如下 012340,30,45,60,90 xxxxxsinx二阶差商一阶差商三阶差商kx()kf x四阶差商34569000057 711866 310.000000.00000.00.0.00000166713811 64470.00.00.0 0 00.006351 713620.00000.000 0 00.000 7264850.0000000.000000 2670.00000000例:已知函数已知函数 的数值表的数值表 sinx00030ix0045iy0.500000.707110600900.866031试求试求 的四次的四次Newton插值多项式，并计算插值多项式

11、，并计算 的近似值。的近似值。sinx0sin40解：解：取节点取节点 ,则则 的差商表如下的差商表如下 012340,30,45,60,90 xxxxxsinx由由Newton插值公式，依次可得插值公式，依次可得1()0 0.01667(0)N xx 21()()0.0000635(0)(30)N xN xxx32()()0.000000726(0)(30)(45)N xN xxxx43()()0.00000000267(0)(30)(45)(60)N xN xxxxx例:已知函数已知函数 的数值表的数值表 sinx00030ix0045iy0.500000.707110600900.866

12、031试求试求 的四次的四次Newton插值多项式，并计算插值多项式，并计算 的近似值。的近似值。sinx0sin40解：解：取节点取节点 ,则则 的差商表如下的差商表如下 012340,30,45,60,90 xxxxxsinx由由Newton插值公式，依次可得插值公式，依次可得1()0 0.01667(4040 0)N 21()()0.0000635(0)(30)N xN xxx32()()0.000000726(0)(30)(45)N xN xxxx43()()+0.00000000267(0)(30)(45)(60)N xN xxxxx0.6668例:已知函数已知函数 的数值表的数值表

13、 sinx00030ix0045iy0.500000.707110600900.866031试求试求 的四次的四次Newton插值多项式，并计算插值多项式，并计算 的近似值。的近似值。sinx0sin40解：解：取节点取节点 ,则则 的差商表如下的差商表如下 012340,30,45,60,90 xxxxxsinx由由Newton插值公式，依次可得插值公式，依次可得1()0 0.01667(4040 0)N 21()()0.0000635(0)4040404(30)00.6414NN32()()0.000000726(0)(30)(45)N xN xxxx43()()+0.0000000026

14、7(0)(30)(45)(60)N xN xxxxx0.6668例:已知函数已知函数 的数值表的数值表 sinx00030ix0045iy0.500000.707110600900.866031试求试求 的四次的四次Newton插值多项式，并计算插值多项式，并计算 的近似值。的近似值。sinx0sin40解：解：取节点取节点 ,则则 的差商表如下的差商表如下 012340,30,45,60,90 xxxxxsinx由由Newton插值公式，依次可得插值公式，依次可得1()0 0.01667(4040 0)N 21()()0.0000635(0)4040404(30)00.6414NN32404

15、040()()0.000000726(0)(404030)(4=5)0.642852NN43()()+0.00000000267(0)(30)(45)(60)N xN xxxxx0.6668例:已知函数已知函数 的数值表的数值表 sinx00030ix0045iy0.500000.707110600900.866031试求试求 的四次的四次Newton插值多项式，并计算插值多项式，并计算 的近似值。的近似值。sinx0sin40解：解：取节点取节点 ,则则 的差商表如下的差商表如下 012340,30,45,60,90 xxxxxsinx由由Newton插值公式，依次可得插值公式，依次可得1(

16、)0 0.01667(4040 0)N 21()()0.0000635(0)4040404(30)00.6414NN32404040()()0.000000726(0)(404030)(4=5)0.642852NN43()()+0.00000000267(0)4040404040(30)(45)(60)40=0.6429588NN0.6668例:已知函数已知函数 的数值表的数值表 sinx00030ix0045iy0.500000.707110600900.866031试求试求 的四次的四次Newton插值多项式，并计算插值多项式，并计算 的近似值。的近似值。sinx0sin40解：解：取节点

17、取节点 ,则则 的差商表如下的差商表如下 012340,30,45,60,90 xxxxxsinx由由Newton插值公式，依次可得插值公式，依次可得04sin(40)(40)=0.6429588N545sin(40 0)(40 30)(40 45)(40 60)(4(40)0 90)5!180R（）5140 10 5 20 505!180 0.000027516 1855,11)(2xxxg例例M3.m55,11)(2xxxg例例M3.m55,11)(2xxxg例例M3.m55,11)(2xxxg例例M3.m55,11)(2xxxg例例M3.m缺点：当插值点个数取得较多时，可能会使得插值多项

18、式产生很大的扰动，进而导致插值余项不收敛，拉格朗日多项式插值的这种振荡现象叫 Runge现象现象55,11)(2xxxg 采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例例M3.m分分段段线线性性插插值值xjxj-1xj+1x0 xnxoynkybyaxfxxx01()已知:函数在节点上的函数值 LI x()若存在折线满足:I xa b(),则称是在上的分段线性插值函数。I xa b(1)(),在上连续；kkI xykn(2)()(0,1,2,),；LkkI xxx1(3)(),在每个小区间上是线性函数。分分段

19、段线线性性插插值值jjjjjjjjjjjxxxxxxxxxlxxxxxx111111,(),0,其 它 xjxj-1xj+1x0 xnxoynjjjIxy lx0()()0()lx1()lx1()ilx()il x1()ilx()nlx0 x1x2x1ixix1ix1nxnx线性无关：显然；1,()0,性质：插值性：节点上的插值基函数：jiiji jjixxl xxx分段线性插值的余项分段线性插值的余项f xC a b2(),定理:设函数，f x()则的分段线性插值余项R xMh21()8 kka x bk nMfxhxx101max(),max()其中 kkkkkkxxxI xlx ylx

20、y111,()()()时明，证：，当故余项kkkfR xf xI xxxxx 1()()()()()()2!kkkkMMhhxx 221,(,)88 kkkhxx1其中。hR x0()0当时，hI xf x0lim()()因而，I xa bf x(),()故在上一致收敛于。66,11)(2xxxg例例用分段线性插值法求插值用分段线性插值法求插值.在在-6,6中平均选取中平均选取41个点作插值个点作插值M4.m比分段线性插值更光滑比分段线性插值更光滑xyxi-1 xiab 在数学上，光滑程度的定量描述是：函数在数学上，光滑程度的定量描述是：函数(曲线曲线)的的k阶阶导数存在且连续，则称导数存在且

21、连续，则称该曲线具有该曲线具有k阶光滑性阶光滑性。三次样条插值三次样条插值三次样条函数满足：三次样条函数满足：在每个子区间在每个子区间 上是一个三次多项式。上是一个三次多项式。(1)S x()iix x 1,(2)S x()在每个内节点上具有二阶连续导数。在每个内节点上具有二阶连续导数。(3)iiS xy()in0,1,.LiiiiiiiS xa xb xc xdxx321(),.三次样条插值三次样条插值S x()S xxxx101(),Sxxx x212(),MnnnSxxxx1(),.例例21(),551g xxx用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值M5.m

22、用用MATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method,extrap)插值方法插值方法被插值点被插值点已知数据已知数据xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ：线性插值；线性插值；spline ：三次样条插值；三次样条插值；cubic ：立方插值。立方插值。pchip ：分段三次分段三次HermiteHermite插值插值缺省时：缺省时：分段线性插值。分段线性插值。注意：所有插值方法都要求注意：所有插值方法都要求x是单调的是单调的例例 在在1-121-12的的1111小时内，每隔

23、小时内，每隔1 1小时测量一次温度，小时测量一次温度，测得的温度依次为：测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/10小时的温度小时的温度值。值。M1.mxy机翼下轮廓线例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。M2.ml 设有一个年产设有一个年产240吨的食品加工厂吨的食品加工厂,需要统计其生产费需要统计其生产费用用,但由于该厂各项资料不全但由于该厂各项资料不全,无法统计。在这种情况无法统计。在这种情况下下,

24、统计部门收集了设备、生产能力和该厂大致相同的统计部门收集了设备、生产能力和该厂大致相同的五个食品加工厂的产量与生产费用资料如下表：五个食品加工厂的产量与生产费用资料如下表：例例M6.m 如何根据已知五个食品加工厂的资料较准确地估计该如何根据已知五个食品加工厂的资料较准确地估计该厂的生产费用呢？厂的生产费用呢？作业6：1、本节内容小结；2、课后作业第3、7题；3、数值实验题第1、3题。二维插值二维插值12004000,12003600.xy二维插值二维插值第一种（网格节点）第一种（网格节点）第二种（散乱节点）第二种（散乱节点）xyO O yx0 0二维插值的定义二维插值的定义已知已知n个节点个节

25、点),.,2,1(),(nizyxiii),(iiyx互不相同，互不相同，构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点,再用再用),(yxf计算插值，即计算插值，即).,(*yxfz),1,0(,),(nizyxfiii即即 注意注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值最邻近插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O 二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。网格节点网格节点将四个插值点（矩形的四个顶点）处的函数值依次简记为：分片线性插值分片线性插值xy (xi,y

26、j)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)O Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4.网格节点网格节点f1f2f3f4分片线性插值分片线性插值xy (xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)O O网格节点网格节点f1f2f3f411()jjijiiyyyxxyxx插值函数为：jiiijjyxxxxyyy)(11)()()(23121jiyyffxxfffx,yf分两片分两片第一片(下三角形区域)：(x,y)满足()f x,yAxByC分片线性插值分片线性插值xy (xi,

27、yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)O O网格节点网格节点f1f2f3f411()jjijiiyyyxxyxx插值函数为：jiiijjyxxxxyyy)(11分两片分两片第二片(下三角形区域)：(x,y)满足)()(),(43141ijxxffyyfffyxf()f x,yAxByC注意注意：(x,y)应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的。双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形如：(,)()()f x yaxb cyd四个待定系数双线性插值双线性插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O

28、网格节点网格节点四个顶点（插值节点）的函数值四个方程四个方程要求：要求：1.x0,y01.x0,y0单调；单调；z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLAB作作网格节点网格节点数据的插值数据的插值插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值2.x2.x，y y可取可取为矩阵，或为矩阵，或x x取行向量，取行向量，y y取为列向量；取为列向量；3.x,y3.x,y的值分别不能超的值分别不能超出出 x

29、0,y0 x0,y0的范围。的范围。例：测得平板表面例：测得平板表面3 3*5 5网格点处的温度分别为：网格点处的温度分别为：82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 82 85 86 试作出平板表面的温度分布曲面试作出平板表面的温度分布曲面z=f(x,y)z=f(x,y)的图形。的图形。1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图先在三维坐标画出原始数据，画出粗糙的温度分布曲图.M7.m2以平滑数据以平滑数据,在在x、y方向上每隔方向上每隔0.2个单位的地方进行插值

30、个单位的地方进行插值.画出插值后的温度分布曲面图画出插值后的温度分布曲面图.通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。M8.m 12004000,12003600.xy第二种（散乱节点）：第二种（散乱节点）：yx0 0 cz=griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作作散点数据散点数据的插值计算的插值计算 要求要求:cxcx取取行行向量，向量，cycy取为取为列列向量向量。被插值点插值方法插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值（双线性插值（缺省时缺省时）cubiccubic 双三次插值双三次插值