数值分析12(插值方法1)

1、1/50*插值的应用背景插值的应用背景拉格朗日插值公式拉格朗日插值公式牛顿插值公式牛顿插值公式插值误差余项插值误差余项Runge反例反例数值分析 122/50*3/50*4/50*趣例趣例1: 图像放大图像放大5/50*趣例趣例2: 工业设计工业设计http:/ Paul de Casteljau 和和Pierre Bzier, 随后美国通用汽车的随后美国通用汽车的其它人一起推动了现在称为三次样条和其它人一起推动了现在称为三次样条和Bzier 样条的建立。样条是通过样条的建立。样条是通过很少的控很少的控制点就能够生成复杂平滑曲线制点就能够生成复杂平滑曲线的方法。的方法。6/50*趣例趣例3:

2、数据可视化数据可视化http:/ 游戏与电影游戏与电影Ref: http:/ 0 1 1)x,y = ginput;n = length(x);s = (1:n);t = (1:.05:n);u = splinetx(s,x,t);v = splinetx(s,y,t);clf resetplot(x,y,.,u,v,-);9/50*数据和插值函数数据和插值函数 如果一个函数如果一个函数P(x)满足满足yi = P(xi) (i=0, n), 那么那么函数函数P(x) 插值了一系列数据点插值了一系列数据点(x0,y0 ), (xn,yn )，其中其中P(x)称为称为插值函数插值函数,点点x0

3、, ,xn称为称为插值节点插值节点。x0 x1x2x3x4xP(x) 函数是描述自然界客观规律的重要工具。函数是描述自然界客观规律的重要工具。10/50*插值问题研究包括如下三个方面插值问题研究包括如下三个方面: :插值函数的选择和构造插值函数的选择和构造插值函数的存在唯一性插值函数的存在唯一性插值误差估计的问题插值误差估计的问题11/50*)()(001010 xxxxyyyxL 过两点直线方程过两点直线方程已知函数表求满足已知函数表求满足: L(x0)=y0, L(x1)=y1的线性函数的线性函数 L(x) x x0 x1 f(x) y0 y1例例 求求 的近似值的近似值(函数值函数值:

4、10.7238)11510.7143)100115(100121101110115 真实值真实值: 10.723812线性插值函数线性插值函数x0 x1(x0 ，y0)(x1，y1)L1(x)可见可见 是过是过 和和 两点的直线。两点的直线。13抛物插值函数抛物插值函数x0 x1x2因过三点的二次曲线为抛物线因过三点的二次曲线为抛物线,故称为抛物插值。故称为抛物插值。 14/50* 选择多项式函数的理由选择多项式函数的理由: 理论方面理论方面多项式函数简单明了的数学性质。多项式函数简单明了的数学性质。有一个简单的原理可以说明什么时候存在给定有一个简单的原理可以说明什么时候存在给定次数的插值多项

5、式。次数的插值多项式。 更重要的是更重要的是计算方面计算方面多项式函数是计算机多项式函数是计算机最基本的函数最基本的函数, 计算多项式函数的值只需用加计算多项式函数的值只需用加和乘运算和乘运算, 且积分和微分均非常方便。且积分和微分均非常方便。插值函数的选择插值函数的选择:15/50*则称则称 P(x) 为为 插值多项式插值多项式, 称称 x0, x1, , xn为为 插插值节点。值节点。 如果如果 P(x)=a0 + a1x + anxn满足满足 P(xk)= yk (k = 0,1,n)考虑区间考虑区间a , b上上(n+1)个点个点a x0 x1xnb。代数插值问题代数插值问题插值条件插

6、值条件16/50*点点,则满足插值条件则满足插值条件 L(xk)= yk (k = 0,1,n)的的 次次数小于等于数小于等于n次次的插值多项式的插值多项式 L(x)=a0 + a1x + anxn存在而且唯一存在而且唯一。 nnnnnnnnnyxaxaayxaxaayxaxaa101111000010证明证明: 由插值条件由插值条件L(x0)= y0L(x1)=y1L(xn)=yn定理定理5.1 若插值结点若插值结点x0, x1, xn 是是 (n+1)个互异个互异17/50*00111 1det( )1nnnnnxxxxAxx 回顾回顾1:非齐次方程组有唯一解的充分必要条件是系数非齐次方程

7、组有唯一解的充分必要条件是系数矩阵可逆矩阵可逆(矩阵可逆的充分必要条件是行列式不等于零矩阵可逆的充分必要条件是行列式不等于零)。系数矩阵行列式不等于零系数矩阵行列式不等于零, 则方程组有唯一解。因此则方程组有唯一解。因此插值多项式插值多项式 L(x) 存在且唯一。存在且唯一。0()0ijij nxx 回顾回顾2: 范德蒙范德蒙(Vandermonde)矩阵矩阵0011102002131101211 1det()()()()()()() ()nnnnnijnnj i nnnxxxxAxxAxxxxxxxxxxxxxxxx 则则( ( ) )= =21()()nnnnxxxx 18/50*注释注释

8、: 插值多项式的存在唯一性说明满足插值插值多项式的存在唯一性说明满足插值条件的多项式存在而且与条件的多项式存在而且与构造方法无关构造方法无关。 只要插值节点互异只要插值节点互异, 则则Vandermonde矩阵矩阵总是非奇异。然而范得蒙矩阵条件数通常很总是非奇异。然而范得蒙矩阵条件数通常很大大,故直接求解方程组是危险的故直接求解方程组是危险的。Hilbert和和Vandermonde条件数条件数:for i=1:10c(i)=cond(hilb(i),2);%vander(1:i)endplot(1:10,c)19/50*)()(001010 xxxxyyyxL 过两点直线方程过两点直线方程已

9、知函数表求满足已知函数表求满足: L(x0)=y0 和和 L(x1)=y1的线性函数的线性函数 L(x)。 x x0 x1 f(x) y0 y120/50*10010101,( )( )xxxxxxxlxxlx 记记x x0 x1l0(x) 1 0l1(x) 0 11100)()()(yxlyxlxL 01010110( )xxxxL xyyxxxx 对称形式对称形式001011110,0,yyyyyy 21/50*二次插值问题二次插值问题 x x0 x1 x2 f(x) y0 y1 y2已知函数表已知函数表求函数求函数 L(x)=a0 + a1x + a2 x2 满足满足: L(x0)=y0

10、 , L(x1)=y1, L(x2)=y2L(x)=l0(x)y0+l1(x)y1+l2(x)y2l0(x) 100l1(x) 010l2(x) 0 01L(x) y0y1y2 xx0 x1 x222/50*)()()(2010210 xxxxxxxxxl 二次插值函数二次插值函数: L(x)=l0(x)y0+l1(x)y1+l2(x)y2，l0(x) 100l1(x) 010l2(x) 0 01L(x) y0y1y2 xx0 x1 x2)()()(2101201xxxxxxxxxl )()()(1202102xxxxxxxxxl 23/50*二次插值基函数图形二次插值基函数图形取取 x0 =

11、0, x1 = 0.5, x2 = 1l0(x)=2(x 0.5)(x 1);l1(x)= 4 x(x 1);l2(x) = 2(x 0.5)x24/50*拉格朗日插值公式拉格朗日插值公式插值条件插值条件:L(xk)= yk (k = 0,1,n)nnnyxlyxlyxlxL)()()()(1100 )()()()()()()(110110nkkkkkknkkkxxxxxxxxxxxxxxxxxl 其中第其中第k 个个插值基函数插值基函数 (k=0,1,，n)。0()()()njkjkjjkxxlxxx 或或0()( )()njnkkjkjj kxxLxyxx 紧紧凑凑格格式式25/50*例例

12、1求插值于点求插值于点(0,2),(1,1),(2,0),(3,-1)的的次数次数小于等于小于等于3的插值多项式的拉格朗日形式。的插值多项式的拉格朗日形式。(1)(2)(3)(0)(2)(3)( )21(01)(02)(03)(10)(12)(13)(0)(1)(3)(0)(1)(2) +01(20)(21)(23)(30)(31)(32) =+2xxxxxxL xxxxxxxx 0()( )()njnkkjkjj kxxLxyxx 26/50*例例2求插值于点求插值于点(-2,-56),(-1,-16),(0,-2),(1,-2),(3,4)的次数小于等于的次数小于等于4的插值多项式拉格朗日

13、形式。的插值多项式拉格朗日形式。(1)(0)(1)(3)(2)(0)(1)(3)( )5616( 21)( 20)( 21)( 23)( 12)( 10)( 11)( 13)(2)(1)(1)(3)(2)(1)(0)(3) 22(02)(01)(01)(03)(12)(11)(10)(13) xxxxxxxxL xxxxxxxxx 23(2)(1)(0)(1) 4(32)(31)(30)(31) =2572xxxxxxx 27/50*程序片段程序片段1:Matlab Code : 多项式插值多项式插值(拉格朗日形式拉格朗日形式)function v = polyinterp(x,y,u)%PO

14、LYINTERP Polynomial interpolation.% v = POLYINTERP(x,y,u) computes v(j) = P(u(j) where P is the% polynomial of degree d = length(x)-1 with P(x(i) = y(i).% Use Lagrangian representation.% Evaluate at all elements of u simultaneously.n = length(x);v = zeros(size(u);for k = 1:n w = ones(size(u); for j

15、= 1:k-1 k+1:n w = (u-x(j)./(x(k)-x(j).*w; end v = v + w*y(k);end28/50*Demo2x=0:3;y=-5 -6 -1 16;u=-.25:.01:3.25;v=polyinterp(x,y,u);plot(x,y,o,u,v,-)symx=sym(x),L=polyinterp(x,y,symx);L=simplify(L);29/50* 拉格朗日形式结构紧凑拉格朗日形式结构紧凑且且形式对称。然形式对称。然而很少用它来计算而很少用它来计算, 这是因为等价的牛顿形这是因为等价的牛顿形式更具操作性而且计算复杂度更低。式更具操作性而且

16、计算复杂度更低。000111 () () ()nnnxyf xxyf xxyf x 牛顿形式相对简单牛顿形式相对简单,但是某些记号首先需但是某些记号首先需要掌握。把数据点看作由某个函数给出要掌握。把数据点看作由某个函数给出, 并并把数据点列成下表把数据点列成下表: :插值多项式的牛顿形式插值多项式的牛顿形式30/50*给定给定x0, x1和和 x2,求二次函数求二次函数 P(x)=a0 + a1(x x0) + a2 (x x0)(x x1)满足条件满足条件 P(x0)=f(x0), P(x1)=f(x1), P(x2)=f(x2) )()()()()()(21202202101011000

17、xfxxxxaxxaaxfxxaaxfa 满足插值条件的关于满足插值条件的关于a0, a1和和a2方程方程插值多项式牛顿形式插值多项式牛顿形式31/50*010101()(),f xf xf xxxx 121212()(),f xf xf x xxx 011201202,f xxf xxf xxxxx 解下三角方程组过程中引入符号解下三角方程组过程中引入符号a0 = f(x0), a1 = fx0, x1, a2 = fx0, x1, x2P(x)= f(x0) + fx1, x2(x x0) + fx0, x1, x2(x x0)(x x1)牛顿插值公式牛顿插值公式:32/50*定义定义5.

18、3 若已知函数若已知函数 f(x) 在点在点 x0,x1,xn 处的处的值值 f(x0), f(x1), , f(xn)。如果如果 i j ,则则( j = 0,1,n-1 ) 一阶差商一阶差商n阶差商阶差商111000,nnnnfxxxxxxxxxff 111()(),jjjjjjf xf xf xxxx 二阶差商二阶差商222111,jjjjjjjjjxxxxxffxxf xx 三阶差商三阶差商123331212,jjjjjjjjjjjjfxxf xxxxxxxxxfx 差商差商(divided difference)33/50*更加一般地考虑更加一般地考虑000110101001 ()(

19、) ()()()()()nnnnnnaf xaa xxf xaa xxaxxxxf x 牛顿插值形式牛顿插值形式求解该方程组可得待定系数如下求解该方程组可得待定系数如下:a0 = f(x0), a1 = fx0, x1, a2 = fx0, x1, x2, ,an = fx0, x1, ,xn。00100101201011( )()()() +()(),)(,),nnf xf xxf xx xfP xxxxxxxxxxxxxxxx 001010121( )()()()+()()()nnP xxxxxxxxxxxaaaxax 34/50*例例3 已知数据如下表已知数据如下表, 试计算数据的差商。

20、试计算数据的差商。 2 561 16 0 2 1 23 4 4014 0 3 13 7 1 2 20111000,nnnnfxxxxxxxxxff 35/50*例例2求插值于点求插值于点(-2,-56),(-1,-16),(0,-2),(1,-2),(3,4)的次数小于等于的次数小于等于4的插值多项式拉格朗日形式。的插值多项式拉格朗日形式。(1)(0)(1)(3)(2)(0)(1)(3)( )5616( 21)( 20)( 21)( 23)( 12)( 10)( 11)( 13)(2)(1)(1)(3)(2)(1)(0)(3) 22(02)(01)(01)(03)(12)(11)(10)(13

21、) xxxxxxxxL xxxxxxxxx 23(2)(1)(0)(1) 4(32)(31)(30)(31) =2572xxxxxxx 36/50*例例4 已知数据如下表已知数据如下表, 试计算数据的插值多项式。试计算数据的插值多项式。 2 561 16 0 2 1 23 4 4014 0 3 13 7 1 2 20 ( )56(2)(2)( +1)(2)( +1)( -0)(2)( +1)( -0)2( -1)14003+P xxxxxxxxxxx (1)(132( -0) ( )5(2 )(40) 6Pxxxx 00100101201011( )()()() +()(),)(,),nnf

22、xf xxf xx xfP xxxxxxxxxxxxxxxx 37/50*例例5 已知数据如下表已知数据如下表, 试计算数据的插值多项式。试计算数据的插值多项式。 2 561 16 0 2 1 23 4 4 7 40 14 0 3 3 13 7 1 0 2 2-1/4 0-9/20 加入一个新的点到加入一个新的点到Lagrange形式所需要额外工作形式所需要额外工作与牛顿形式进行比较是很有趣的。牛顿形式具有与牛顿形式进行比较是很有趣的。牛顿形式具有Lagrange形式所缺少的形式所缺少的”实时更新实时更新”性质。性质。38/50*011011(),()()()()nknkkkkkkknf xf

23、 xxxxxxxxxxx 差商的性质差商的性质00100101201011( )()()() +()(),)(,),nnf xf xxf xx xfP xxxxxxxxxxxxxxxx 0()( )()()njkkjkjj kxxP xf xxx 差商的值不依赖于差商的值不依赖于x0,x1,xn的次序。的次序。39/50*压缩的概念压缩的概念: 观测的离散数据可以想象成现实中无穷多观测的离散数据可以想象成现实中无穷多信息的代表。通过给定数据求出插值函数意味信息的代表。通过给定数据求出插值函数意味着用简单的规则代替无穷多信息。尽管期待这着用简单的规则代替无穷多信息。尽管期待这种简单规则精确地反映

24、实际情况是不现实的种简单规则精确地反映实际情况是不现实的, 但是它但是它可以充分接近实际可以充分接近实际。 这一类压缩是这一类压缩是有损的压缩有损的压缩, 即它会产生误差。即它会产生误差。用简单规则代替无穷多信息时会产生多大的误用简单规则代替无穷多信息时会产生多大的误差差, 这是我们下面研究的内容。这是我们下面研究的内容。40/50*00111010)(yxxxxyxxxxxL 两点线性插值两点线性插值插值误差余项插值误差余项: R(x) = f(x) L1(x) 由插值条件知由插值条件知 R(x)=C(x) (x x0)(x x1)即即 f(x) L(x) = C(x) (x x0)(x x

25、1) C(x) = ?41/50*Rolle过山车过山车:( ),( , ),( )( ),( )=0f xa ba bf af bf 设设函函数数在在闭闭区区间间连连续续 在在开开区区间间可可微微且且则则至至少少存存在在一一点点 使使得得。42/50*回顾回顾: 拉格朗日中值定理拉格朗日中值定理( ),( , ),( , )( )( ) ( )=f xa ba ba bf bf afba 设设函函数数在在闭闭区区间间 连连续续 在在开开区区间间可可微微则则至至少少存存在在一一点点使使得得 。( )( ):( )( )( )()f bf aF xf xf axaba 思思路路43/50*Ln(

26、x)是满足是满足Ln(xk)= f(xk) 的的n次插值多项式次插值多项式, 则则对任何对任何xa, b, 在在(a, b)内存在一点内存在一点 使得使得(1)1( ( )( )( )(1)!nnnnfxRxf xLxnx 余余项项101( )()()()nnxxxxxxx 。其中其中定理定理5.2 设设 f (x) 在在a, b连续且在连续且在(a, b) 具有具有n+1阶导数阶导数, x0,x1,xn是是a, b内互不相同的节点内互不相同的节点。( ) x 44/50*证明证明: 记记 n+1(x) =(x x0) (x xn)(0, ),()(),( , )kkkxxknf xL xa

27、b 首首先先注注意意到到如如果果则则在在内内任任意意选选择择 均均成成立立。(0, ),kxxkna bt 如如果果则则考考虑虑在在 内内定定义义关关于于 的的函函数数如如下下。00()()( )( )( ) ( )( )()()nntxtxg tf tL tf xL xxxxx 00()()()()()() ( )( )()()() =0 ( )( )00kkkkknkkkknt xxxxxxxg xf xL xf xL xxxxxxxf xL x 对对 = =有有45/50*00()()( )( )( ) ( )( )()() = ( )( ) ( )( )0nnnnnnt xxxxxg

28、xf xLxf xLxxxxxf xLxf xLx 对对 = = 有有(1)( )2Rolle,( )( , )1, ( )( , )1ng tng ta bngta b 所所以以至至少少有有 + + 个个互互异异的的零零点点。根根据据定定理理 在在内内至至少少存存在在 + + 个个相相异异的的零零点点。以以此此类类推推在在内内至至少少存存在在 个个零零点点。(1)(1)0(1)0( ),( +1)!0( )( ) ( )( )()()( )( )( )()()( +1)!nnnnnnnxngff xLxxxxxff xLxxxxxn 故故存存在在使使得得整整理理即即得得46/50*(1)0(1)10( )()Taylor( )=Taylor()( +1)!( )()( +1)!)=nnnnR xfxxxxnfxnR xx 插插值值余余项项与与展展式式余余项项的的异异同同插插值值余余项项展展式式余余项项注释注释: 47/50*例例1 给出如下数据给出如下数据 x 0.32 0.34 0.36 sin(x) 0