数值分析 matlab方法 插值法

1、第五章第五章理学院应用数学系理学院应用数学系立体化教学资源系列立体化教学资源系列数值分析数值分析数值分析数值分析插值法插值法5.1 5.1 引言引言人口普查数据人口普查数据(千人千人) 年194019401950195019601960197019701980198019901990人口132,165132,165151,326151,326179,323179,323203,302203,302226,542226,542249,633249,633推测推测19301930年、年、19651965年、年、20102010年的人口年的人口. . 194019451950195519601965

2、197019751980198519901.21.41.61.822.22.42.6x 105 美国人口预测美国人口预测数值分析数值分析插值法插值法一天一天2424小时零点开始每间隔小时零点开始每间隔2 2小时环境温度（度）小时环境温度（度）: :12 9 9 10 18 24 28 27 25 20 18 15 1312 9 9 10 18 24 28 27 25 20 18 15 13， 温度预测温度预测 推测中午推测中午1 1点（即点（即1313点）的温度点）的温度. .05101520255101520253002040608010012014016018020002468101214

3、1618数值分析数值分析插值法插值法数据如表数据如表: 机翼上缘轮廓曲线机翼上缘轮廓曲线 x 0.00 4.74 9.50 19.00 38.00 57.00 76.00y 0.00 5.32 8.10 11.97 16.15 17.10 16.34x95.00 114.0 133.0 152.0 171.0 190.0y14.63 12.16 9.69 7.03 3.99 0.00数值分析数值分析插值法插值法020406080100120140160180200024681012141618三次样条函数画出机翼曲线三次样条函数画出机翼曲线(MATLAB ) x=0.00 4.74 9.50

4、19.00 38.00 57.00 76.00 95.0 114.0 133.0 152.0 171.0 190.0; y=0.00 5.32 8.10 11.97 16.15 17.10 16.34 14.63 12.16 9.69 7.03 3.99 0.00; xx=0.0:0.1:190; yy=interp1(x,y,xx,spline); plot(x,y,*), pause, hold on, plot(xx,yy)数值分析数值分析插值法插值法利用利用MATLAB函数函数peaks产生一个山顶曲面数据产生一个山顶曲面数据 山顶曲面山顶曲面 x,y,z = peaks(10);me

5、sh(x,y,z),hold on,plot3(x,y,z,r*),hold off -4-2024-4-2024-10-50510数值分析数值分析插值法插值法通过插值作出更加精细的山顶曲面通过插值作出更加精细的山顶曲面 figure(2)xi,yi = meshgrid(-3:.1:3,-3:.1:3);zi = interp2(x,y,z,xi,yi); mesh(xi,yi,zi) -4-2024-4-2024-10-50510数值分析数值分析插值法插值法美国大学生数学建模竞赛美国大学生数学建模竞赛1986年年A题题. 水道测量问题水道测量问题 水道测量数据表水道测量数据表/m x 12

6、9.0 140.0 108.5 88.0 185.5 195.0 105.5 y 7.5 141.5 28.0 147.0 22.5 137.5 85.5 z 1.22 2.44 1.83 2.44 1.83 2.44 2.44 x 157.5 107.5 77.0 81.0 162.0 162.0 117.5y -6.5 -81.5 3.0 56.5 -66.5 84.0 -38.5z 2.74 2.74 2.44 2.44 2.74 1.22 2.74数值分析数值分析插值法插值法水道地貌形状图水道地貌形状图 水道深度等值线图水道深度等值线图 50100150200-1000100200-3

7、-2.5-2-1.5-1-0.5XYDepth-2.9-2.8-2.8-2.8-2.7-2.7-2.7-2.6-2.6-2.6-2.5-2.5-2.5-2.5-2.5-2.4-2.4-2.4-2.4-2.4-2.3-2.3-2.3-2.3-2.3-2.2-2.2-2.2-2.2-2.2-2.2-2.1-2.1-2.1-2.1-2.1-2.1-2-2-2-2-2-2-1.9-1.9-1.9-1.9-1.9-1.8-1.8-1.8-1.8-1.8-1.7-1.7-1.7-1.7-1.7-1.7-1.6-1.6-1.6-1.6-1.6-1.5-1.5-1.5-1.5-1.5-1.4-1.4-1.4-

8、1.4-1.3-1.3-1.3-1.3-1.2-1.2-1.2-1.1-1.1XY80100120140160180-50050100MATLAB程序见下页程序见下页 数值分析数值分析插值法插值法MATLAB程序程序（数据表中水深：（数据表中水深：z值加负号）值加负号）： x=129.0 140.0 108.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5;y=7.5 141.5 28.0 147.0 22.5 137.5 85.5 -6.5 -81.5 3.0 56.5 -66.5 84.0 -38.5;z=-1

9、.22 2.44 1.83 2.44 1.83 2.44 2.44 2.74 2.74 2.44 2.44 2.74 1.22 2.74;xi=linspace(min(x),max(x),40);yi=linspace(min(y),max(y),40);Xi,Yi=meshgrid(xi,yi);Xi,Yi,Zi=griddata(x,y,z,Xi,Yi,cubic);mesh(Xi,Yi,Zi);xlabel(X);ylabel(Y);zlabel(Depth);figure(2);c,h=contour(Xi,Yi,Zi,-1.1:-0.1:-2.9);clabel(c,h);xlab

10、el(X);ylabel(Y);求求?)245sin( ),1 , 0(nixi)(xpniyxpii,1 , 0,)( 正弦函数表如下正弦函数表如下 x643xysin212223本章问题本章问题数值分析数值分析插值法插值法0 xnxx1x0yxysin1yny函数值计算函数值计算设设 , ,且已知且已知 上上 个点个点)(xfy,ba , a b1n),1 , 0(niyi插值问题插值问题插值区间插值区间被插值函数被插值函数 )(xfy 插值函数插值函数 )(xp),1 ,0(nixi插值节点插值节点 插值条件插值条件 )(xp为多项式为多项式.插值多项式插值多项式: 已知数据表，求已知数

11、据表，求f(x)的近似函数的近似函数.求简单函数求简单函数 ,使满足使满足的对应函数值为的对应函数值为niyxpii,1 , 0,)(a,b存在且唯一存在且唯一. . 所谓所谓插值问题插值问题：就是已知：就是已知 被插值函数在插值区间上一些被插值函数在插值区间上一些互异节点的函数值互异节点的函数值, ,求插值函求插值函数数 , , 使 满 足 插 值 条 件使 满 足 插 值 条 件)(xp满足满足 个互异节点条件个互异节点条件 1n), 1 , 0()(niyxpiinnnnxcxcxccxp2210)(数值分析数值分析插值法插值法的多项式的多项式【定理【定理1 1】的求法的求法: : )(

12、xp1)1)待定系数法待定系数法( (解方程组解方程组) )2)2)构造法构造法: :LagrangeLagrange、NewtonNewton两点对称式方程：两点对称式方程：点斜式方程：点斜式方程：几何意义：几何意义：过两个已知点过两个已知点，求直线方程，求直线方程 1 , 0),(iyxii数值分析数值分析插值法插值法)()(0010101xxxxyyyxL101001011)(yxxxxyxxxxxL01011010)(,)(xxxxxlxxxxxl一、线性插值一、线性插值(n=1) 由两点式由两点式 ：11001)()()(yxlyxlxL插值基函数：插值基函数： 1 , 0, 0,

13、1)(ikikikxlkiik5.2 5.2 拉格朗日插值拉格朗日插值5.2.1 线性插值与线性插值与抛物插值抛物插值212021012101200201021)()()()()()(yxxxxxxxxyxxxxxxxxyxxxxxxxx数值分析数值分析插值法插值法(n=2) 二、抛物插值二、抛物插值几何意义：几何意义： 过三个点过三个点,2, 1 ,0),(iyxii求抛物线求抛物线 基函数性质基函数性质 2211002)()()()(yxlyxlyxlxLnkkknyxlxL0)()(nk, 1 ,0nkiiikikxxxxxl0)()(2 , 1 , 0, 0, 1)(ikikikxlk

14、iik5.2.2 n次拉格朗次拉格朗日插值日插值数值分析数值分析插值法插值法例例1 1 已知函数已知函数表表?245sin解解1 1）线性插值线性插值 60 x41x取两点取两点，则，则 2264621464)(1xxxL0.603553)245(245sin1L抛物插值抛物插值 23)43)(63()4)(6(22)34)(64()3)(6(21)36)(46()3)(4()(2xxxxxxxL0.609577)247(247sin2Lx643xysin2122232 2）0.60876140.6087614 )245sin(抛物比线性插值精确抛物比线性插值精确. . 【注】【注】【注】【注

15、】数值分析数值分析插值法插值法 计算上不方便计算上不方便. .计算机算法容易实现计算机算法容易实现. )()!1()()()()()1(xwnfxLxfxRxnnn【定理【定理2 2】 ,baxnjjxxxw0)()(其中其中, ,)()!1()(1xwnMxRnn（1 1）误差估计）误差估计 （2 2） 节点的位置、个数节点的位置、个数 有关有关 n余项与余项与 1nMx、（3 3）当）当 是是 的多项式时的多项式时 )(xfn)()(xfxLn这是因这是因为为0)()1(xfn)(max)1(),(1xfMnbaxn5.2.3 5.2.3 插值余项与误差估计插值余项与误差估计缺点：缺点：L

16、agrange插值优点：插值优点： 公式整齐对称，适合理论推导公式整齐对称，适合理论推导,时，有时，有1)(xfnkkxl01)(数值分析数值分析插值法插值法例例2 2 估计例估计例1 1中中 与与 的误差的误差. . )245 (1L)245(2Lxxfsin)()sin()(xxf )cos()(xxf )(! 2)(1021xxxxMxR)(max),(210 xfMxxx 4,6x0.0061)4245)(6245(! 24sin)245(1R)()(! 3)(21032xxxxxxMxR)(max),(320 xfMxxx 3,6x0.00097)3245)(4245)(6245(!

17、 36cos)245(2R解解，有，有1 1） 线性插值的误差估计线性插值的误差估计. .因为因为其中，其中，所以，所以2 2） 抛物插值误差估计抛物插值误差估计. .因为因为其中，其中，所以，所以由由function yi=Lagrange(x,y,xi)%Lagrange%Lagrange插值多项式，其中，插值多项式，其中，% %x x为向量，全部的插值节点；为向量，全部的插值节点；% %y y为向量，插值节点处的函数值为向量，插值节点处的函数值; ;% %xixi为标量或向量，被估计函数的为标量或向量，被估计函数的 自变量；自变量；% %yiyi为为xixi处的函数估计值处的函数估计值.

18、 .n=length(x); m=length(y);% %插值点与它的函数值应有相同个插值点与它的函数值应有相同个 数数. .if n=merror(The lengths of X and Y must be equal!);return;endLagrangeLagrange插值的插值的MATLABMATLAB程序：程序：LagrangeLagrange.m.myi=zeros(size(xi);for k=1:nw=ones(size(xi);for j=1:k-1 k+1:n% %输入的插值节点必须互异输入的插值节点必须互异If abs(x(k)-x(j)epserror(the D

19、ATA is error!);return;Endw=(xi-x(j)./(x(k)-x(j).*w;endyi=yi+w*y(k);end数值分析数值分析插值法插值法调用调用插值插值函数函数Lagrange.mLagrange.mx=0:3;y=-5 -6 -1 16;x=0:3;y=-5 -6 -1 16;xi=-.25:.01:3.25;xi=-.25:.01:3.25;yi=yi= Lagrange Lagrange (x,y,1.5) (x,y,1.5)yi=yi= Lagrange Lagrange (x,y,xi); (x,y,xi); plot(x,y,o,xi,yi)plot

20、(x,y,o,xi,yi)( )5, 6, 1,16.f x 求求f(x)的的LagrangeLagrange插值多项式，插值多项式，程序计算得，程序计算得， yi = -4.6250. yi = -4.6250. -0.500.511.522.533.5-10-50510152025已知当已知当x=0,1,2,3时时, ,数值分析数值分析插值法插值法例例3并求并求f(1.5)(1.5)的近似值的近似值.ikikiiikiiikiiixxxxxfxxxfxxxf,11211【定义定义1 1】已知互异节点数据已知互异节点数据 ), 1 , 0)(,(nixfxii一阶差商一阶差商 iiiiiix

21、xxfxfxxf111)()(,二阶差商二阶差商 iiiiiiiiixxxxfxxfxxxf212121,K 阶差商阶差商 差商性质差商性质 (1)(1)K 阶差商可由节点函数值的线性组合表示阶差商可由节点函数值的线性组合表示,10kxxxfkjkjjjjjjjjxxxxxxxxxxxf01110)()()()(5.3.1 5.3.1 差商的概念和性质差商的概念和性质5.3 5.3 差商与牛顿插值差商与牛顿插值数值分析数值分析插值法插值法（2 2）差商的对称性：对调节点，差商值不变）差商的对称性：对调节点，差商值不变. . )!1()(,)1(10nfxxxxfnn（3 3） 例例4 4 己知

22、函数表，求四阶差商己知函数表，求四阶差商. . 0 06 61 110104 43 34646181814/3 14/3 4 4828236366 61/31/36 6212212 656529/329/311/1511/151/151/15ix)(ixf,1iixxf,3iixxf,4iixxf,21iiixxxf数值分析数值分析插值法插值法,1010kijkjixxxxxfxxxxxf设所求多项式为设所求多项式为 )()()()()(110102010nnnxxxxxxaxxxxaxxaaxN根据根据n+1+1个插值条件来确定系数个插值条件来确定系数 naaa,100 xx00ya增加一个

23、节点增加一个节点 nx, ,只需由只需由 加上一个加上一个n次项次项)(1xNn1xx ,1001011xxfxxyya22)(yxNn11)(yxNn2xx ,2102xxxfa)(xNn得到得到插值多项式为插值多项式为 )()()( )()()()()(1101110102010nnnnnnxxxxxxaxNxxxxxxaxxxxaxxaaxN5.3.2 5.3.2 牛顿插值公式牛顿插值公式数值分析数值分析插值法插值法【注注】ix)(ixf0 x1x2x3x)(0 xf)(1xf,10 xxf)(2xf,21xxf,210 xxxf)(3xf,32xxf,321xxxf,3210 xxxx

24、f 一阶差商一阶差商 二阶差商二阶差商三阶差商三阶差商步骤步骤1 1：构造差商表：构造差商表NewtonNewton插值方法插值方法 步骤步骤2 2：牛顿插值公式：牛顿插值公式).()(,)(,)(,)()(110210102100100nnnxxxxxxxxxxfxxxxxxxfxxxxfxfxN数值分析数值分析插值法插值法 1/2=0.5 1/2=0.5 =0.707107 =0.707107 0.7910900.791090 =0.866025 =0.866025 0.607024 -0.3515390.607024 -0.351539例例5 5 x643xysin2122231 1）

25、构造差商表构造差商表 ix)(ixf6432223二阶二阶一阶一阶2 2）线性插值）线性插值 603553. 0)245(245sin1N抛物插值抛物插值 )245(245sin2N609577. 0)4245)(6245(351539. 0)245(1N数值分析数值分析插值法插值法解解)()!1()()()()()1(xwnfxLxfxRxnnn)(,)()()(10 xwxxxxfxNxfxRnnn)!1()(,)1(10nfxxxxfxnn这是因为这是因为 【注注】对于数表函数对于数表函数f(x),差商形式的余项公式差商形式的余项公式 更具一般性更具一般性. .5.3.3 5.3.3 插

26、值余项插值余项与误差估计与误差估计数值分析数值分析插值法插值法牛顿插值的牛顿插值的MATLABMATLAB程序：程序：Newton_interpNewton_interp.m.mfunctionfunction yi,Y=Newton_interp(x,y,xi)yi,Y=Newton_interp(x,y,xi)%Newton%Newton插值多项式，其中，插值多项式，其中，% %x x为向量，全部的插值节点为向量，全部的插值节点; ;% %y y为向量，插值节点处的函数值为向量，插值节点处的函数值; ;% %xixi为标量，被估计函数的自变量；为标量，被估计函数的自变量；% %yiyi为为

27、xixi处的函数估计值；处的函数估计值；%Y为差商表为差商表.n=length(x); m=length(y);%插值点与它的函数插值点与它的函数值应有相同个数值应有相同个数.if n=merror(The lengths of X and Y must be equal!);return;end数值分析数值分析插值法插值法% %计算差商表计算差商表. .Y=zeros(n);Y(:,1)=y;Y=zeros(n);Y(:,1)=y;for k=1:n-1for k=1:n-1 for i=1:n-k for i=1:n-k% %输入的插值节点必须互异输入的插值节点必须互异. . if abs

28、(x(i+k)-x(i)epsif abs(x(i+k)-x(i)eps error(the DATA is error(the DATA is error!);error!); return; return; end endY(i,k+1)=(Y(i+1,k)-Y(i,k+1)=(Y(i+1,k)-Y(i,k)/(x(i+k)-x(i);Y(i,k)/(x(i+k)-x(i); end endendend % %计算计算NewtonNewton插值公式插值公式N(xi).N(xi).yi=0;yi=0;for i=1:nfor i=1:n z=1; z=1; for k=1:i-1 for

29、k=1:i-1 z=z z=z* *(xi-x(k);(xi-x(k); end end yi=yi+Y(1,i) yi=yi+Y(1,i)* *z;z;end end 数值分析数值分析插值法插值法等距节点等距节点 ), 1 , 0(0nkkhxxk步长步长h 已知已知 )(kkxfy 【定义【定义2 2】一阶向前差分一阶向前差分 kkkyyy11kkkyyy)2()2(2121hxfhxfyyykkkkk一阶向后差分一阶向后差分 中心差分中心差分 二阶向前差分二阶向前差分 kkkyyy12m阶差分阶差分 kmkmkmyyy111111kmkmkmyyy211211kmkmkmyyy5.4 5

30、.4 差分与等距节点插值差分与等距节点插值 5.4.1 5.4.1 差分的概念和性质差分的概念和性质数值分析数值分析插值法插值法差分的基本性质差分的基本性质 差商与差分的关系差商与差分的关系 00101101,yhxxyyxxf02221021,yhxxxf010!1,yhkxxxfkkk差分与导数关系差分与导数关系 kkkhkxxxfxf!,)(100kkkkhfhkkf)(!)()()()()()(0kkkfhxf),(0kxxx位于数据表的表头位于数据表的表头 差商用差分代替差商用差分代替 ，则，则NewtonNewton公式简化为等距节点插值公式简化为等距节点插值 5.4.2 5.4.

31、2 等距节点插值公式等距节点插值公式数值分析数值分析插值法插值法（1）（2）一、牛顿向前差分插值公式一、牛顿向前差分插值公式 等距节点等距节点 ), 1 , 0(0nkkhxxk) 10(0tthxxhktxxk)( hktxxk)( )()(,11010kkxxxxxxxxxf)1() 1()0(!10kththtyhkkk)1() 1(!10ktttykk牛顿差商插值多项式牛顿差商插值多项式 )()(,)(,)(,)()(110210102100100nnnxxxxxxxxxxfxxxxxxxfxxxxfxfxN牛顿向前差分插值多项式牛顿向前差分插值多项式 ) 1() 1(!2/ ) 1(

32、)(002000ntttnyttytyythxNnn余项余项 ),(),() 1()!1()()(01)1(nnnnxxnttthnfxR数值分析数值分析插值法插值法x位于数据表的表尾位于数据表的表尾 牛顿插值公式牛顿插值公式 )()(,)(,)(,)()(11011211xxxxxxxxxfxxxxxxxfxxxxfxfxNnnnnnnnnnnnnnn)01(tthxxn牛顿向后差分插值公式牛顿向后差分插值公式 ) 1() 1(!) 1(! 2)(2ntttnyttytyythxNnnnnnnn余项余项 ),(),() 1()!1()()(01)1(nnnnxxnttthnfxR例例7 xy

33、cos) 1 . 0, 4 , 1 , 0(hkkhxk求求cos0.048cos0.048及及cos0.35cos0.35的的近似值并估计误差近似值并估计误差.1） 首先构造差分表首先构造差分表数值分析数值分析插值法插值法解解二、牛顿向后差分插值公式二、牛顿向后差分插值公式 0.00.00.10.10.20.20.30.30.40.41.000001.000000.99500 -0.005000.99500 -0.005000.98007 -0.01493 -0.009930.98007 -0.01493 -0.009930.95534 -0.02473 -0.00980 0.00013 0

34、.95534 -0.02473 -0.00980 0.00013 0.92106 -0.03428 -0.00955 0.00025 0.000120.92106 -0.03428 -0.00955 0.00025 0.00012ixiy)( yy )(22yy )(33yy )(44yy 2 2）用牛顿前插公式计算）用牛顿前插公式计算cos0.048 cos0.048 048. 0 x1 . 0h48. 0/ )0(hxt 2/ ) 148. 0(48. 000993. 048. 0005. 01)048. 0()048. 0cos(4 N6/ )248. 0)(148. 0(48. 000

35、013. 099884. 024/ )348. 0)(248. 0)(148. 0(48. 000012. 07-554101.0921)4)(3)(2)(1(! 5)048. 0(ttttthMR0.3894)4 . 0sin(5M误差估计为误差估计为其中其中 数值分析数值分析插值法插值法function function yi,Y=Newton_forward(x,yi,Y=Newton_forward(x,y,xi)y,xi)%Newton%Newton前插公式，其中，前插公式，其中，%x%x为行向量，等距插值节点；为行向量，等距插值节点；%y%y为行向量，插值节点处的为行向量，插值节点

36、处的函数值函数值; ;%xi%xi为标量，被估计函数的为标量，被估计函数的自变量；自变量；%yi%yi为为xixi处的函数估计值；处的函数估计值；%Y%Y为差分表为差分表. . for k=1:i-1for k=1:i-1 z=z z=z* *(t-k);(t-k); end end% %计算初始值计算初始值. .h=x(2)-x(1);t=(xi-x(1)/h;h=x(2)-x(1);t=(xi-x(1)/h;% %计算差分表计算差分表Y.Y.n=length(y);n=length(y);Y=zeros(n);Y(:,1)=y;Y=zeros(n);Y(:,1)=y;for k=1:n-1

37、for k=1:n-1Y(:,k+1)=diff(y,k);Y(:,k+1)=diff(y,k);zeros(k,1);zeros(k,1);endend% %计算计算NewtonNewton前插公式前插公式N(xi).N(xi).yi=Y(1,1);yi=Y(1,1);for i=1:n-1for i=1:n-1 z=t; z=t; yi=yi+Y(1,i+1)yi=yi+Y(1,i+1)* *z/prod(1:i);z/prod(1:i);endend牛顿向前差分插值的牛顿向前差分插值的MATLABMATLAB程序：程序：Newton_forward.mNewton_forward.m数值

38、分析数值分析插值法插值法function yi,Y=Newton_backward(x,y,xi)function yi,Y=Newton_backward(x,y,xi)%Newton%Newton后插公式，其中，后插公式，其中，%x%x为行向量，等距插值节点；为行向量，等距插值节点；%y%y为行向量，插值节点处的函数值为行向量，插值节点处的函数值; ;%xi%xi为标量，被估计函数的自变量；为标量，被估计函数的自变量；%yi%yi为为xixi处的函数估计值；处的函数估计值；%Y%Y为差分表为差分表. .% %计算初始值计算初始值. .n=length(x); h=x(n)-x(n-1);t

39、=(x(n)-xi)/h;n=length(x); h=x(n)-x(n-1);t=(x(n)-xi)/h;% %计算差分表计算差分表Y.Y.牛顿向后差分插值的牛顿向后差分插值的MATLAB程序：程序：Newton_backward.m数值分析数值分析插值法插值法n=length(y); Y=zeros(n);Y(:,1)=y;n=length(y); Y=zeros(n);Y(:,1)=y;for k=1:n-1for k=1:n-1 Y(:,k+1)=zeros(k,1);diff(y,k); Y(:,k+1)=zeros(k,1);diff(y,k);endend% %计算计算Newto

40、nNewton后插公式后插公式N(xi).N(xi).yi=Y(n,1);yi=Y(n,1);for i=1:n-1for i=1:n-1 z=t; z=t; for k=1:i-1 for k=1:i-1 z=z z=z* *(t-k);(t-k); end end yi=yi+Y(n,i+1) yi=yi+Y(n,i+1)* *(-1)i(-1)i* *z/prod(1:i);z/prod(1:i);endend数值分析数值分析插值法插值法对于对于例例7，调用插值函数，调用插值函数Newton_forward.m求插值求插值.程序计算得，程序计算得，cos0.048的近似值为的近似值为yi

41、 = 0.9988 x=0:0.1:0.4;y=1 0.995 0.98007 0.95534 0.92106;xi=0.048;yi,Y=Newton_forward(x,y,xi)再调用插值函数再调用插值函数Newton_backward.m求插值求插值. x=0:0.1:0.4;y=1 0.995 0.98007 0.95534 0.92106;xi=0.35;yi,Y=Newton_backward(x,y,xi)程序计算得，程序计算得，cos0.35的近似值为的近似值为yi = 0.9394 数值分析数值分析插值法插值法例如，设计飞机外形曲线是密合程度更好的插值例如，设计飞机外形曲线

42、是密合程度更好的插值. .HermiteHermite插值问题插值问题 已知数据表已知数据表 ix0 x1xnx)(iixfy 0y1yny)(iixfm0m1mnm求插值多项式求插值多项式)(xH，满足，满足 nimxHyxHiiii, 1 , 0,)(,)(由由2n+1个条件，可唯一确定一个次数不超过个条件，可唯一确定一个次数不超过2n+2的多的多项式项式 1212221012)(nnnxcxcxccxH5.5 5.5 埃尔米特插值埃尔米特插值数值分析数值分析插值法插值法假设有两个多项式假设有两个多项式 及及反证法反证法: : 唯一性唯一性)(12xHn)(12xHn )()()(1212

43、xHxHxQnn次数不超过次数不超过2 2n+1+1且满足且满足 ., 1 , 0, 0)(, 0)(nixQxQii), 1 , 0(nixxi都是都是Q( (x)的二重零点的二重零点 共有共有2 2n+2+2个零点个零点 0)(xQ)()(1212xHxHnn)(),(xxkk), 1 , 0(nk基函数基函数为待定的为待定的2n+1次多项式次多项式 nkkkkknmxyxxH012)()()(存在性存在性的构造法证明的构造法证明满足满足5.5.1 5.5.1 埃尔米特插值多项式的存在唯一性埃尔米特插值多项式的存在唯一性数值分析数值分析插值法插值法0)(, 0, 1)(ikkiikxiki

44、kxikikxxkiikik, 0, 1)(, 0)(nik, 1 , 0,， )()()(),()()(222211xlbxaxxlbxaxkkkk0)(2)(, 1)(111kkkkkkkxlaxbxax)(21kkxlakkkkxxlxab)(21111)()()(21 )(2xlxlxxxkkkkk)()()(2xlxxxkkkHermite插值多项式插值多项式 nkkkkkkkknxlmxxyxlxxxH0212)()()()(21)(数值分析数值分析插值法插值法210)22(12)()()!22()()(nnnxxxxxxnfxR),(ba其中其中且与且与x有关有关 1n11001

45、1003)()()()()(mxmxyxyxxH2010011120110100)(21 ()(,)(21 ()(xxxxxxxxxxxxxxxxxx201011201100)()(,)()(xxxxxxxxxxxxxx),(,)()(! 4)()(102120)4(3xxxxxxfxR插值余项插值余项 5.5.3 5.5.3 三次埃尔三次埃尔米特插值多项式米特插值多项式5.5.2 5.5.2 埃尔米特埃尔米特插值余项插值余项数值分析数值分析插值法插值法,)1(21 )(,) 1)(21 ()(2120 xxxxxx,) 1()(,) 1()(2120 xxxxxx1100113)()()()

46、(mxmxyxxH5403. 0) 1() 1(8415. 0)23(222xxxxxx47821. 0)5 . 0(3H2） 估计误差估计误差 2120)4(3)()(! 4)()(xxxxfxR8415. 01sinsinmax)(max)4(4xxfM0022. 0) 15 . 0()05 . 0(! 48415. 0)5 . 0(223R数值分析数值分析插值法插值法例例8 8 )5 . 0(3H，并估计误差，并估计误差. 已知已知y= =f( (x) )的函数值及导数值如表，求次数不超的函数值及导数值如表，求次数不超 )(3xP1133)(),2 , 1 , 0( ,)(fxPiyxP

47、ii，满足，满足1）求插值多项式）求插值多项式 )(3xP通过三点的二次插值通过三点的二次插值 )(,)(,)(102100100 xxxxxxxfxxxxfxf)(,)(,)()(1021001003xxxxxxxfxxxxfxfxP)()(210 xxxxxxA 113)(fxP由条件由条件)()(,)(210101210101xxxxxxxxxfxxfxfA2） 求插值余项求插值余项 )()(! 4)()()()(2210)4(3xxxxxxfxPxfxR),(20 xx其中其中 且依赖于且依赖于x. 数值分析数值分析插值法插值法例例9过过3的多项式的多项式解解来确定来确定A。数值分析数

48、值分析插值法插值法5.6.1 5.6.1 高次插值高次插值的病态性质的病态性质5.6 5.6 分段低次分段低次插值插值龙格龙格(Runge) 5 , 5)1 (1)(2xxf10n)10, 1 , 0(5kkxk)1 (12xy-5-4-3-2-1012345-0.500.511.52- 1/(1+x2)- L10(x)-5-4-3-2-101234500.10.20.30.40.50.60.70.80.91原 函 数分 段 线 性 插 值高次插值多项式不稳定高次插值多项式不稳定 节点节点)(10 xLy 在在-5,5上的图形上的图形 通常不用高次插值，而用分段低次插值。通常不用高次插值，而用

49、分段低次插值。 【注注】二、分段三次二、分段三次Hermite插值插值一、分段线性插值一、分段线性插值1111)(kkkkkkkkhyxxxxyxxxxxI12121)(21 ()()(21 ()(kkkkkkhyhxxhxxyhxxhxxxI12121)()(kkkkkkmhxxxxmhxxxx,1kkxxxkkxxh11, 1 , 0nk其中其中其中其中,1kkxxx1, 1 , 0nk一、分段线性插值一、分段线性插值 MhxIxfxRh8)()()(2)(max,max)4(41xfMxxhbxakkk其中其中数值分析数值分析插值法插值法5.6.2 5.6.2 分段低次插值方法分段低次插

50、值方法5.6.3 5.6.3 分段低次插值余项分段低次插值余项)(max,max1xfMxxhbxakkk 二、分段三次二、分段三次Hermite插值插值 44384)()()(MhxIxfxRh其中其中)1 (1)(2xxf)10, 1 , 0( ,5kkxk给定函数给定函数，取插值节点，取插值节点 分段三次分段三次HermiteHermite插值（右图）插值（右图）分段线性插值（左图）分段线性插值（左图）-5-4-3-2-101234500.10.20.30.40.50.60.70.80.91-5-4-3-2-101234500.10.20.30.40.50.60.70.80.91数值分析

51、数值分析插值法插值法bxxxan10:)(kkxfy ), 1 , 0(nk1） 插值条件：插值条件： 2） 分段条件：分段条件： 3） 光滑条件：光滑条件： nkyxSkk, 1 , 0,)(,1kkxx1, 1 , 0nk332210)(xaxaxaaxS,)(2baCxS则称则称S(x)三次样条插值函数三次样条插值函数. 【定义【定义2 2】5.7 5.7 三次样条插值三次样条插值5.7.1 5.7.1 三次样条插值三次样条插值已知已知求函数求函数S(x)满足满足: 确定确定S(x)需补充两个边界条件需补充两个边界条件.数值分析数值分析插值法插值法（1） 第一种边界条件第一种边界条件nn

52、fxSfxS)(,)(00nnfxSfxS )(,)(00自然边界自然边界 （2） 第二种边界条件第二种边界条件（3） 第三种边界条件第三种边界条件0)()(0 nxSxS nmmxSxS02 , 1 , 0m 例例10己知己知1) 1 (, 0)0(, 1) 1(fff2n.1 , 0,)(,0 , 1,)()(11213110020300 xdxcxbxaxsxdxcxbxaxsxS1) 1 (, 0)0(, 0)0(, 1) 1(1100ssss 求求f(x)在自然边界条件下的三次样条插值多项式在自然边界条件下的三次样条插值多项式.解解 待定系数待定系数 . 1, 0, 0, 11111

53、0000cbaddcba数值分析数值分析插值法插值法)0()0(),0()0(1010ssss 1010,bbcc0) 1 (, 0) 1(10 ss026 , 0261100baba 联立联立 0,23,2110101010ddccbbaa.1 , 0,2321,0 , 1,2321)(2323xxxxxxxSxkkMS（x）在节点在节点处的处的弯矩弯矩待定待定而构造出而构造出S(x)的表达式的表达式kkMxS )(三弯矩法就是由三弯矩法就是由) 1, 1 , 0(nk,1kkxx)(xS 数值分析数值分析插值法插值法5.7.2 5.7.2 三弯矩法三弯矩法是一次式是一次式,)(111 kk

54、kkkkkkxxxhxxMhxxMxS线性插值线性插值数值分析数值分析插值法插值法kkkxxh1其中其中 积分两次，得到积分两次，得到 112212)(2)()(CMhxxMhxxxSkkkkkk2113316)(6)()(CxCMhxxMhxxxSkkkkkk利用插值条件利用插值条件 11)(,)(kkkkyxSyxS可求出常数可求出常数21,CC )(62)(2)()(111221kkkkkkkkkkkkMMhhyyMhxxMhxxxSkkkkkkkkkkkhxxhMyMhxxMhxxxS121331)6(6)(6)()(. 1, 1 , 0,)6(1211nkxxxhxxhMykkkkk

55、kkkM)0()0(kkxSxS确定确定，利用连续条件，利用连续条件数值分析数值分析插值法插值法kkkkkkkkkkkkkkhyyMhMhhyyMhMh111111163631, 2 , 1211nkdMMMkkkkkk,6,11111kkkkkkkkkkkkxxxfdhhhhhh（1 1）对第一种边界条件）对第一种边界条件 nnfxSfxS)(,)(00),(6),(621020212111101001011211nnnnnnnnnxxffhddfxxfhMMM三对角方程组三对角方程组 其中其中求弯矩的求弯矩的 数值分析数值分析插值法插值法nnnnnnnnnMdddMdMMMM1122011

56、122112222120222 的三对角方程组的三对角方程组 （2 2）对第二种边界条件）对第二种边界条件 nnfxSfxS )(,)(00121,nMMM求求 （3）对第三种边界条件）对第三种边界条件 nnnnnnnnnMdddMdMMMM1122011122112222120222121,nMMM求求 三对角三对角方程组方程组数值分析数值分析插值法插值法【注注】)(xS)(xS)(xS ,ba)(xf)(xf )(xf 及及于于上均分别一致收敛于上均分别一致收敛于其中其中)(xS)(xS)(xS 收敛最快，收敛最快，次之，次之，最慢。最慢。 例例1111 已知函数，已知函数， 求三次样条插

57、值函数求三次样条插值函数 )(xS满足边界条件满足边界条件 14)2(,75. 0)5 . 1(SS对于第一种或第二种边界条件对于第一种或第二种边界条件 2 , 1 , 0,)(max)()(max4)4()()(ihxfCxSxfibxaiiibxankxxhhhkkkknk, 1 , 0,max1108/3,24/1,384/5210CCCx-1.5 0 1 2-1.5 0 1 2y0.125 -1 1 90.125 -1 1 95.7.3 5.7.3 三次样条插值的误差估计与收敛性三次样条插值的误差估计与收敛性kkkd,kM1 1）根据给定数据和边界条件计算）根据给定数据和边界条件计算，

58、写出确定，写出确定的线性方程组的线性方程组 36186 . 6621005 . 025 . 0004 . 026 . 000123210MMMM2 2）解方程组，得）解方程组，得 16, 4, 4, 53210MMMM3 3）将所得）将所得kM代入代入各子区间上各子区间上, , 所求三次样条插值函数所求三次样条插值函数).21 (3642),10(12),05 . 1(12)(23223xxxxxxxxxxS解解一、一、interp1函数函数命令格式：命令格式：yi=interp1(x,y,xi,method) 1、linear分段线性插值分段线性插值命令格式：命令格式：yi=interp1(

59、x,y,xi,linear)画画y=sin(x)y=sin(x)在在0,100,10曲线曲线, ,取曲线上插值节点取曲线上插值节点例例12kxk10, 1 , 0k折线图折线图.MATLAB.MATLAB计算如下计算如下及其函数值，画出分段线性插值及其函数值，画出分段线性插值x=0:10;y=sin(x);xi=0:.25:10;yi=interp1(x,y,xi);x=0:10;y=sin(x);xi=0:.25:10;yi=interp1(x,y,xi);plot(x,y,o,xi,yi,k-,xi,sin(xi),k:)plot(x,y,o,xi,yi,k-,xi,sin(xi),k:)，xlabel(x),ylabel(y)xlabel(x),ylabel(y)5.8.1 5.8.1 一维插值函数一